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INTRODUCTION 

The trend in propeller aircraft has been toward increasing the speed, 
size, and' horsepower. As a result, there is an increasing demand for the 
design of propellers which are efficient and yet produce minimum noise. 

This requires accurate determinations of both the flow over the propeller 
blade surfaces and the acoustic field induced by the moving propeller. 

Much effort has been devoted recently to the development of a more 
sophisticated propeller theory. This theory has proceeded from the early 
simple momentum model of W. J. M. Rankine (1, 1865) and R. E. Froude (2, 

1889) and the blade-element model of W. Froude (3, 1878) and S. Drzewiecki 
(4, 1920) to the vortex theory first proposed by F. W. Lanchester (5, 1907), 
the lifting-line model of S. Goldstein (6, 1929) and, finally, to the lifting 
surface model of H. Ludwieg and I. Ginzel (7, 1944). One of many important 
advancements in lifting-line theory is Lerbs' calculation of the induction 
factors (8, 1952), which allows the velocities at each blade section to be 
determined with great accuracy. This important calculation, plus other more 
sophisticated mathematical models, e.g., P. C. Pien (9, 1961), J. E. Kerwin 
(10, 1964), and W. B. Morgan (11, 1968), makes it possible today to design 
a propeller based on fluid dynamic principles. 

One of the important problems of aeroacoustics is the determination 
of the sound from a rotating propeller. Historically, L. Gutin (12, 1936) 
was the first to theoretically investigate this sound for a static rotating 
propeller, using an equivalent distribution of dipoles in the propeller 
disk. His method later was extended and generalized by I. E. Garrick and 
C. E. Watkins (13, 1953) to the case of an in-flight propeller by considering 
the pressure dipoles that represent the thrust and torque force to be sub- 
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jected to a uniform rectilinear motion. The general theory of aerodynamic 
sound given by M. J. Lighthill. (14, 1952 and 15, 1954) has been extended 
to situations containing arbitrarily moving boundaries in unbounded space 
by J. E. Ffowcs Williams and D. L. Hawkings (16, 1969), using the theory 
of generalized functions. The surface is replaced by a discontinuity in 
the flow-field, around which the motion of the fluid medium is assumed to 
be known. Other important works in rotational propeller noise are given 
in Refs. 17-25. In addition, K. Karamcheti and Y. H. Yu (26, 1974) 
have studied the hovering rotor propeller, minimizing the far field inten- 
sity subject to aerodynamic constraints. 

This paper is concerned with the design of propellers for minimum 
noise. The paper is divided into three parts. In order to relate aero- 
dynamic propeller design and propeller acoustics, the first part includes 
the necessary approximations and assumptions involved, the coordinate sys- 
tems and their transformations, the geometry of the propeller blade, and 
the - problem formulations including the induced velocity, required in 
the determination of mean lines of blade sections, and the optimization 
of propeller noise. The second part is devoted to the numerical formula- 
tion for the lifting-line model. The third part presents some applications 


and numerical results. 
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PART 1.. BASIC THEORY J\ND FORMULATION ' 

A. Assumptions and Consequences 

An exact propeller design is not possible in any theoretical analysis 
so that a number of simplifying assumptions must be made. Except in certain 
special cases like stall-flutter, where nonlinearity is of essence, the 
aerodynamic tools should be mathematically linear to ensure the possibility 
of finding a solution with reasonable time and effort. The following as- 
sumptions, based on this concept, are generally made in theoretical pro- 
peller design and in the calculation of the sound pressure due to the 
moving propeller blades. 

-1. The propeller is operating in an unlimited stationary fluid’ with 
a constant advance velocity and a constant angular velocity. 

2. The fluid is inviscid. Although all real fluids are compressible 
to a greater or lesser extent, under normal conditions the effects of com- 
pressibility are unimportant at low speed, and consequently“the density of 
the fluid will be assumed to be constant in developing the vortex theory. 
However, from the acoustic viewpoint, the compressibility of fluid is im- 
portant so that the fluid is restored to be compressible in the acoustic 
formulation. 

3. The propeller consists of a set of identical, symmetrically spaced 
blades attached to a hub. The hub effect is ignored so that it is not 
necessary to satisfy the hub boundary conditions. 

4. The blade sections are thin and the blade is not heavily loaded. 

In this case the disturbance velocities produced by the propeller are small 
compared with the propeller advance velocity and rotational velocity. There- 
fore, the deviation between the blade surfaces and the stream surfaces formed 
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by the relative undisturbed flow is small. This assumption permits us 
to treat the problem as a logical extension of linearized finite wing 
theory where the tangency condition is satisfied on the mean line of the 
profile. It also enables separation of the loading and thickness effects. 
However, from the acoustic viewpoint, this idea implies that the quad- 
ruple sources, the Lighthill stress, are negligible, since they contain 
only those second-order perturbation terms which are dropped upon line- 
arization. 

5. Each propeller blade is replaced by a reference surface which 

is the projection of the actual blade outline on the helical surface with 
pitch angle the hydrodynamic advance angle obtained from the lifting 
line theory. A distribution of bound vortices for loading effects and 
sources and sinks for thickness effects are placed upon this reference 
surface. The vortices are distributed in both the chordwise and spanwise 
directions. The variation of. strength of vorticity necessitates free 
vortices being shed from the bound vortices. These free vortices form 
helical surfaces behind the propeller and extend to infinity in the pro- 
peller-fixed coordinate system. 

6. Upon neglecting the quadrupole sources, the blade loading and 
thickness are the only acoustic sources that will be considered. 

7. The effects of slipstream contraction and centrifugal force on 
the shape of the free vortex sheets are ignored. Consequently, each of 
the free vortices has a constant diameter and constant pitch downstream 
which may be varied along the radius. 

8.. Body forces are ignored. 

9. The two-dimensional chordwise pressure distributions are preserved 
in the three-dimensional flow. 
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We shall write all the quantities used in the formulation in non- 
dimensional form by referring all velocities to a reference velocity, 

(V may be chosm to be the advance velocity of the propeller) , and by 
referring all linear dimensions to a characteristic length, R^, (pro- 
peller radius) . The pressure and the force per unit area are made non- 
dimensional with respect to V , where 5* 0 is the undisturbed density 
of the fluid and time is referred to 1/H , where il is' the angular velo- 
city of the propeller. Also, the circulation is non -dimens ionali zed with 
respect to 27 fRpVp, and the strengths of the vortex sheets are referred 
to V . Further, the expression 
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is referred to as the reference advance coefficient, which is the advance 

coefficient of the propeller if V is chosen to be the advance velocity 

P 

of the propeller. Dimensional values are denoted by primes, so that, for 
example 



( 2 ) 


B. Coordinate Systems 

Two main coordinate systems, shown in Figs. 1 and .2 , are used 
in the analysis. One is the "space-fixed" coordinate system, which has a 
rectangular (x 1 , x 2 , x 3 ) coordinate system (x-system) , a rectangular (y 1 , 
y 2 , y 3 ) coordinate system (y-system) , and a spherical (§,{S),5) coordinate 
system -system) . The other one is the 'propel ler- fixed" coordinate system. 
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which has a rectangular (x, y, z) coordinate system (z-system) , a cylin- 
drical (x, r, Q ) coordinate system (r-system), and a curvilinear (s, n, r) 
coordinate system (s -system) . Also defined are the unit tangent vectors, 
e w , to the coordinate curves w, where w = x^ y^, r, , and so forth. 

Propeller-Fixed Coordinate System: 

All the propeller-fixed coordinates are attached to the propeller, 
translating and rotating with the propeller. 

• 1. A rectangular (x, y, z) coordinate system (z-system) 

x-axis = axis of revolution of propeller with positive 

distance measured downstream 

y-axis = selected so as to pass through the tip of one 

blade 

z-axis = selected so as to complete the right-handed 

system 

Z = position vector of a space point referred to the 

center of the z-system 

2. A cylindrical (x, r, $ ) coordinate system (r-system) 

x-axis = defined as before 

r-coordinate = radial coordinate 

9 -coordinate = angular coordinate, measured clockwise starting 
from the y-axis when looking downstream 

3. An orthogonal curvilinear (s, n, r) coordinate system (s-system) 
r-coordinate = defined as before 

s-coordinate - a helix whose non-dimensional pitch is 


?. = 2?rr to* fair) - ztt A.i(r) 


(3) 
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where jS^ (h) is the hydrodynamic advance angle, 
and AiCr) is the hydrodynamic advance coefficient. 
Furthermore, the s-coordinate is the intersection 
of the reference surface representing the first 
blade with a circular cylinder of radius r, 
measured downstream along the helix 
n-coordinate = selected so as to complete the right-handed system 

Ipace-Fixed Coordinate System: 

1. A rectangular (x^, x 2 , x^) coordinate system, referred to as an 
observer system (x-system) 

x-j-axis = axis of revolution of the propeller with positive 

distance measured downstream 

The x^-axis, x^-axis, and x^-axds are fixed in space 
and complete a right-handed system. At time t=0, the 
origin of the x-system coincides witlrthat of the 
z-system and the x 2 ~axis and y-axis make an angle Q o 
measured clockwise from the x 2 ~axis when looking 
downstream. 

X = position vector of an observer referred to the 

x-system 

2. A rectangular (y^, y 2 , y^) coordinate system, referred to as a 
source system (y-system) 

- The y^-axis, the y 2 -axis, and the y^-axis are se- 

lected so as to coincide with the x^-axis, the x 2 ~ 
axis, and the x 3 -axis, respectively. 

Y = position vector of an acoustic source referred to 
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the y - system 

3. A spherical coordinate system, also referred to as an 

observer system (^-system) 

-coordinate = distance measured from the origin of the x-system 
to the observer 

$ -coordinate = angular coordinate, measured clockwise starting fron 
x 2 ~axis when looking downstream 
(3) -coordinate = angular coordinate, measured from the x^-axis 

f ' -x 


C. Some Useful Coordinate Transformations 

Following are some useful coordinate transformations which will be 


used in the formulation of the problem. Also given are some relations 
among the unit tangent vectors of the different coordinate systems: 


1. x~ system and ^ -system 

X, = § C4HS, (g) 

Xz — % c<rt 5 (B) 

Xj — § -A**' £ 0 

2. 2 -system and y-system 

it - X - X-f Vp -i: 

$ 2 =■ C*t- 0 a ) -f -g A*** (■£ - do) 

y - C't-Qo') + * ^ C-t- do) 


e * ! 


i o 0 ] 


f e 3i | 

e v 

~ 

o CA<(,t-& 0 ) - 


% 

s 


0 (-t- $») cen(-b-0 o )j 


\%l 


(4) 


(5) 


( 6 ) 



11 


3. z- system and r- system 


% - X 
$ = r c^c 0 

i = r ^ 0 




f, 0 

0 1 




— 

o c*x>0 




e *i 


^ 0 ydiri# 

C-cri- 0j 


, e ‘l 


1 


f 1 0 o 1 


e * 

e r 

— 

o Cd-Z0 ^4i*v 0 


e * 

l e ‘l 


0 zM*'® Q&L& 


e *l 


( 7 ) 


( 8 ) 


4. y- system and r- system 


$1 = “A “* -t 

f Wt( 0 4 00 *t) C9) 

J 3 = r A** (.0 *Q 0 ~±) 


U 


' / 0 

0 ^ 


f i 

! v 

e r 

s 

o c^t C ^ -f- Oo ~ t ) 

>4*»v ( 0 -t 0e> •"■£) 






(0-f- 0o-t) j 


1% 
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5. r-system and s-system 

r = r 

, . Ai.' (r)x + r\»- $*) 

Jr 2 ~ "-^Cr) 

jtr — v A*Jv) £ § * Sk") 

n - — — 

J r 2 + 

where S k — & -coordinate of the point at the tip of the kth blade. For a 
symmetrical blade arrangement, these angles are 


k=i,2,...,$ ( 12 ) 


D. Geometrical Considerations 

Figure 3 shows a projected view of the propeller, looking in the down- 
stream direction (positive x) . The angular coordinates & L and 9- of the 
lea din g and trailing edges, respectively, define the projected blade outline. 

The reference surface representing the kth blade is the projection of 
the actual blade outline on the helical surface: 

H k (x,r 3 e) = x - Ac(r) ( & - S K ) ~o (133 

An expanded view of the. s-n plane showing a typical blade section 
oriented approximately along the s-coordinate curve is illustrated in Fig. 4. 
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y 



Fig. 3 Projected blade outlines; 3-bladed 


propeller. 
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Fig. 4 Illustration of blade section. 


15 


Here s L (r) and s T (r) are the s-coordinates of the leading and trailing 
edges, while the total expanded chord length is l(r). The incidence d (r) 
of the blade section at radius r is defined as the angle between the 
chord line of the section and the s-axis. The incidence is considered 
to be positive when the chord line has greater pitch than the reference 
surface. The mean line c(r,s) and' the thickness t(r,s), are also shown. 

It is noted that the positioning of the actual section in the n-direction 
is immaterial since the blade section will be represented by singularities 
distributed along the s-axis. 

The pitch of each of the reference surfaces is illustrated in Fig. 5. 
For a lightly load ed propeller, the strictly linearized case, these refe- 
rence surfaces will coincide with those swept out by the undisturbed rela- 
tive flow past the radial lines 


6 - Sk 

X - 0 

through the tip of each blade. 

It is seen from Fig. 5 that the non-dimensional pitch is 


(14) 


Pjr) = Zv .^ r = 2tc r tamper) = 2ir A f os) 

where = advance angle of the propeller 

Ap = advance coefficient of the propeller 

Therefore, for the lightly loaded propeller theory, it is sufficient 

to set P. (r) = P (r) . For the moderately loaded propeller, the nonlinear 
l o 

problem being approximated by an equivalent linear one with the considera- 
tion of perturbation velocities (induced velocities), the pitch of the 
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reference surface is increased over the lightly loaded case. From Fig. 5 
it is seen that 



zirVpiVf ajt ul(rt) 
Rf - a t 


£7rr( Vf -t aa.(r)) 

^ + uZw 


= 2 7T f tanji^Y) 
- 2 T 


where 7\^(r) = hydrodynamic advance coefficient 
PiW = hydrodynamic advance angle 

u*(r) = axial component of induced velocity from lifting-line 
theory 

u*(r) = tangential component of induced velocity from lifting 
line theory 

The effective inflow velocity is 

- W + _ Xf» * 

For convenience in the analysis which follows , information about the 
reference surface is given: 

1. #+**) = coordinate of any point on the kth reference 
surface, expressed in the r- system. <p is the angular coordinate 
of the corresponding point on the first reference surface. 

2 . The unit tangent vector at ( (sH-f, $+& has three compo- 


v% 


<£(l) 


nents as follows: 
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_e r = c«.(4+^(c) ^ 

~ Ji^Al(f) { Ai(f)e *“ S >*■(.* 

, CIS) 

e " 'Jf j +A*(f) { $ e x + Mtlux’Lt*^' '<*} -XJ. f)o4(M*5 e^J 

3. 44= J f +' A?(?) dfdf Ciaj 

= infinitesimal area element of the reference surface at 

(Ms)<t>> ?,<!>+ Sk) 

4. ds * J J 1 -V A*(?) df (20) 

= infinitesimal line element along the helix at 

f , 4> + s* ) 



It is noticed that £ } <£> are the dummy cylindrical coordinates of 
the point on the first reference surface. 

H. Mathematical Formulation 

This section is concerned with the formulation of the problem. It 
is divided into three subsections. The first one will be referred to as 
the aerodynamic formulation, dealing with the calculation of the mean 
lines of the blade sections. The second one will be referred to as the 
acoustic formulation, dealing with the acoustic problem caused by the 
moving propeller. The third subsection is concerned with formulation of 
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the criteria for optimization of propeller noise. 

E-l. Aerodynamic Formulation 

The mean line of the blade section relative to the helical chord at 
radius r is determined by the relative induced velocity normal to the 
reference surface (V^^(r,<j>) - V*(r)), where i- s the total . 

induced velocity obtained from lifting-surface theory and V* (r) is the 
induced velocity obtained from lifting-line theory. We shall formulate 
this problem following closely the work of Kerwin and Leopold (10, 1964) , 
based on the assumptions made in Section A. 


Lifting Surface Theory: 

1. Vortex Distribution 

The total bound circulation around the blade section at radius r 
will be defined as P (r) so that 

LO) = J 4j»(r,s) ds = V(t)j hr, s) 

- 2 v V*cr) P(r) 


( 22 ) 


where L(r) 
Ap(r,s) 


= non-dimensional total lift force per unit radius, L'Cr^^Y^R^ 

= non-dimensional pressure differential due to velocity dis- 

2 

continuity, Ap Vj^V 


r(r,s) = non-dimensional strength of the radially oriented bound vor- 
tices that induce a discontinuity in the streamwise velocity 
of + 1/2 r(r,s) at each point on the reference surface, 
2T(r,s)/Vp. 

$ 

P (r) = non-dimensional circulation, P(r)/27fRpVp. 



To satisfy the Helmholtz lav/ of continuity of vortices, this system of 
bound vortices must be accompanied by a system of trailing vortices whose 


axis is in the e direction along a helix of pitch P . (r) . 

s * 

If the strength of the helical vortex sheet is defined as Y s (r) be- 
hind the trailing edge, then 


KM 


= - zir 


dVdr) 

dr 


^ d (23) 


With this expression and Eq. (22)-, 
from the blade are found to be 

s T Cr) 

= *■(*■>*> Js 


the strengths of the free vortices shed 

9f0") 

=-ij' tit, 6)J r 2 +*l(r) 16 
e L Vr) 


-- j ^{H«e>fFTjpT}<l0 

i w — 

— Jr C r, 0 t 1 i r 1 -* AiCr) 


(24) 

lit. 
d r 


The first term on the- right is due to the radial change in the bound 
vortices. The second term is due to the change of ^(r) with respect to 
r along the leading edge. It follows from this equation that, within 
the reference surface,. the free vortex strength £((r, 0) can be expressed 
as follows: 
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It is evident that 


y s ( r,e) = -2vdpH 


07 - ^ (26) 


2. Source Distribution 

The thickness of the blades can be represented by a source sheet 
distribution. The connection between the source strength ^(rjS), 
which induces a discontinuity in normal velocity +_(r( r ,s) at any point 
on the surface, the effective inflow velocity V*(r) and the local change 
of the thickness with chordwise dimension is 


<r(r,s) 


vV> 


3tcns) 

•as 


(27) 


where t(r,s) = non-dimensional thickness as shorn in Fig. 6 — 

0-(r,s) = non-dimensional source strength, (T' , (r,s)/Vp 

3. Induced Velocities 

As mentioned before, it is the main purpose of this section to com- 
pute the induced velocities due to loading and thickness. The computation 
of the induced velocities at points on the reference surfaces representing 
the blades of the propeller enables us to determine the way in which the 
blade sections should be cambered and oriented with respect to the effec- 
tive inflow if a propeller is designed with a prescribed pressure loading 
and thickness. 


3 . 1 Loading 




thickness form 
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From the Biot-Savart law, the induced velocities and at 

any point PCA;.^)^/ V, 8 ) on the first blade by the bound vortices and 
trailing vortices, respectively, are found to be 


8 




<^r ^ 0 


7T D 3 ^ 


(28) 


.1 ,<» 


B 


e £ x D 


y ( 1p)=jj 

?«&*»&< f) 


dA 


(29) 


- Cb ^ -Cfc) 

V(f) - V (£) + V(f) 


(30) 


where r^ = hub radius 

d ={A i (r)a~x^?)4 > }e x 4* ^ it e -i* c ^ -fr W } -t- { e - 5^} e ^ 

the vector distance from source point f, <£ + S(«) 

to the reference point P (A^W&) f, 0 ) on the reference surface. 

D = }d I 

y(t>) _ ^ n( j uce( ^ velocity due to bound vortices 

(P) = induced velocity due to trailing vortices 
V*^ (P) = total induced velocity due to loading 

Upon submitting Eqs. (18) and (19) into Eqs. (28) and (29)., evaluating the 
vector product, and converting velocity components into cylindrical coordinates. 


we have 
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<b> 

u a (r,e) = 


r &r(« g 

♦fJ j 

■Mi. +-9 t (f) 


K^Si* hz£L. j f*+AZ(?) ™ 


.f 0 t(?) 


j XMi£ 




jf 2 + A^(?) (32) 


U r (ne)=“j J KW) 1“ ^ J S’ + Ale f) < 33} 


l CO 




-^j j iti! 


(34) 


UtM^) )W^— s — 3 ^C 3 S) 

f^*=e L (?) 


Lt) 

. u r cce) 


4-7T j 


Ax(f) (f +5 k- 0 ) + (A;Cr) 9 - Ai(?J 4} f e*% (f + Sk-9) 

P 3 


d4JJ( 36) 
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where D =(^18-^)^+ r'+f -2 (37) 

u *- b -* (r, 9 ) = non-dimensional axial component of induced velocity due to 

3 

bound vortices 

u ^ (r, 9 ) = non-dimensional tangential component of induced velocity 

due to bound vortices 

u ^(r,0 ) » non-dimensional radial component of induced velocity due 

- to bound vortices 

u ^ (r, © ) = non-dimensional axial component of induced velocity due 

to trailing vortices 

u t W (r, 0) = non-dimensional tangential component of induced velocity 

due to trailing vortices 

u ^ (r, 6) = non-dimensional radial component of induced velocity due 

to trailing vortices 

* 

It should be noticed that the expressions for the induced velocity 
due to trailing vortices are different from those presented by Kerwin 
by a factor ( j* 2 + 

3.2 Thickness 

The velocity induced at any point P ( AxCD 0; T > & ) on the first 
reference surface by sources distributed over B blades is found by taking 
the gradient of the velocity potential of the sources 

— r l 8 i 

v7P)= j «■(?.*> ^ ™ 


Upon substituting Eq. (19) into Eq. (38), taking the gradient, and 
converting to cylindrical coordinates, one has 



26 







r^3 


where 
f si 

u (r , $ ) = non-dimensional axial component of induced velocity due to 
thickness 

u t ^ s ^(r, 0 ) = non-dimensional tangential component of induced velocity 

due to thickness 
fsl 

u r v (r>£) = non-dimensional radial component of induced velocity due to 
thickness 

.3.3 Induced Velocity Normal to Blade 

As mentioned before, in order to calculate the mean lines of the blade 
sections, the normal component, u , of the total induced velocity due to 
loading and thickness is required at any point P ( A^(\r) 0 ^ j r , 6 } on the 

first blade. This normal component of the induced velocity is related to 
the axial and tangential components of the total induced velocity at that 
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point by the expression: 


. —/kt) — cs) s 

u*(r,e)=(v (ne)+V (ne>)*e* 


= u n (n*) + u„,(r, 0) + UyvCiTj^) 


( 42 ) 


(i) _ct) 

c»;« = v ftaj-e* 


rUg(ng) - A^W^tCr, 6) 
J r" + a*(D 


(43) 


<i> 

u*(na) « 


Ct) aa, 

*_ ru«.(r J e)-A < cCr)UtOG 0) 
V ( T, 0 ) * <- n - : 

J IT* * A^Crt 


(44) 


(cj (s) 

u n Cne) = y (r,«)-e n = 


<s> te) 

rugftfl -Xi® t^e ere) 
J r* •+ a*o-) 


(45) 


4. Mean Line Calculation (Boundary Conditions) 

The mean line of the blade section at radius r can be obtained by 
considering the boundary conditions at the blade surface. The boundary 
condition is that the resultant velocity at any point on the blade surface 
must be tangent to the surface at that point. In the case of linearized 
theory, the surface to which the resultant velocity is tangent, at a given 
radius and chordwise position, is defined to be the surface tangent 
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to the mean line at that point as shown in Fig. 7. 

From Fig. 7 and within the concept of linearized theory , the boundary 
condition at any point on the blade surface is 


3S 


Cm ( n $) 


UnCr, £) - u£(h ) 

v*w 


(46} 


where 

c (r,s) = ordinate of the mean line of the blade section at radius 

r relative to the chord, measured in the direction 
starting from the helical chord. 

u*(r) = normal component of induced velocity from lifting-line 

theory 

Upon integrating Eq. (46) , we have 


■Sf (b) 


Cm cr,s) = - f js 

J |/W 


(47) 


Cm (r, S) 

J r*+ 


& T (n *, t 

j - Un fr? ^ 


(48) 




V T Or) 


Introducing the camber c(r,s) and the ideal angle of attack cij (r) , 
Eq.' - (47) may be expressed as 

5r(h) 

, , , ... ( — Ciftlt) r/icn 

c(ns)+ (s r -s)±a*.c/ A (jr) s ds 

' \/(r) 
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To determine o^fr), this integral is evaluated from the leading edge 
to the trailing edge and c(r,s) is taken to be zero. Having obtained the 
mean line, the blade section is obtained by adding the thickness to it. 

Lifting-Line Theory: 

As described in the preceding section, the total induced velocity • 
normal to the blade surface can be (divided into two parts. One is the 
velocity, u* (r) , due to the lifting -line. helical vortex sheet model. The 
other is the velocity resulting from spreading out the concentrated . 
vortex lines to the desired blade outline and adding thickness to the 
blade. In order to obtain the second part, u*(r) must be obtained first. 

The assumptions underlying the theory of moderately loaded_ propellers 
permit one to design a propeller either to produce a given thrust or to 
absorb a given power, and to compute the ideal thrust and power and, there- 
fore, efficiency, utilizing only the lifting -line representation of the 
propeller . 

1. Lifting-Line Induced Velocity 

As mentioned before, we are concerned with a hubless lifting-line 
propeller so that only the relationship between circulation and lifting 
line disturbance velocity for a hubless propeller is determined. 

Since the lifting- line model of the propeller is a degenerate case 
of the lifting-surface model, the induced velocity components may be 
obtained from Eqs. (31) through (36). The induced velocity due to thick- 
ness is ignored. With the assumptions that the propeller consists of a 
set of identical, symmetrically..spaced blades attached to a hub, the induced 
velocity components associated with thickness and chordwise distribution 
on the blade disappear, leaving only the wake term. Introducing the 
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lifting-line induction factors (8, 1952), we have the induced velocity 
components at the first lifting line 


icx(r) = 


i f' Jrc?) 

zj dS 


r I ±£*J) 
k ?- r 


) <V 


( 50 ) 
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u t (r) 
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(52) 


U$(r) 


ui(ir) <&Lfa(r) - t4^.(Y)+AUfo(r) 

( Jn(nf) v , p 
2 J d f f - r ) a J 




(53) 


V*(r) = tfj'to e x + u*cr) e e + u*(r) 

where 

u*(r) = lifting-line axial component of induced velocity 

u| (r) * lifting-line 'tangential component of induced velocity 

u*(r) = lifting-line radial component of induced velocity 

u*(r) = lifting-line normal component of induced velocity 

I (r, f ) - lifting-line axial induction factor 

a ' 
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) = lifting-line tangential induction factor 
I^(r,£) = lifting-line radial induction factor 
I (x, J) - lifting-line normal induction factor 
^ = Cauchy principal value integral 

The lifting line induction factors are defined as 


ijr ,?), fi .(Wir -tricky*)} d 

J k=i n* d y 
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(54) 
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* J p * 3 
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where 


{f 2 + r 2 - 2 ?r +j*- 2 ll(f) j 


Vz 


(58) 
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/■ 


A »(r) = 


r { v f + uttn] 

r f + 


(59) 


It is evident from Eq. (53) and the expressions for the induction factors 
that the components of induced velocity can not easily be obtained from Eq. 
(53) for a given circulation P(r) since all induction factors are functions 
of X^(r) which, in turn, depends on the axial and tangential components 
of the induced velocity through Eq. (59). In order to obtain these components 
of induced velocity, an_ iterative scheme, which will be discussed in the 
second part, must be applied. A brief summary of the evaluation of induc- 
tion factors may he found in Appendix A. 


2. Thrust and Power Coefficients 

Shown in Fig. 8' are the velocity and force components at each radius 
according to the theory of moderately loaded propellers. The lift L(r) 
can be expressed in terms of bound circulation P (r) (see Eq. 22) as 


L (» = 2 ir l Ar) PW 


(60) 


The effective inflow velocity V*(r) can be written in terms of u* 

a 

u| , A 3 and ^(r) as follows 


V*( r) = 


“V (f) 
M — 


(61) 
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or 


v*00 


■f 

Vp •* UaCr ) 
(rt 


(62) 


From Fig. 8, we obtain 


RLOO= LCrt { 1 - 6<r ^ ’t*+piCri } 


(63) 


^ , 6CW , 

F t cir)= L0r)s4i*fc(r) { t + ‘^^ r) J (64) 

where 

F (x) = axial component of force, per unit radius acting on the 

Ct * 

blade section at r 

F t (r) = tangential component of force, per unit radius acting on 
the blade section at r 

6 (r) = drag/ lift ratio of the blade section at r 

Upon substitution Eqs. (60)’, (61), and (62) into Eqs. (63) and (64), 
and integrating both Eqs. (63) and (64), we have the thrust and power coef- 
ficients, non-dimensionalized with respect to the reference velocity and 
the radius of the propeller, as follows 

C T - , — r -4B frwf-t + 

7 J tA f 1 
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2 

ipv 



-+ 

-tanker) 



( 66 ) 


where 


T = total propeller thrust 
P = total power absorbed 

The efficiency of the propeller is the ratio of power output to the 
power input which in this case is simply 



(67) 


Finally, the ideal thrust coefficient, the ideal power coefficient, 
and the ideal efficiency are obtained directly from Eqs.'(65),(66) and (67) 
by setting £ (r) equal to zero: 



( 68 ) 



(69) 



( 70 ) 
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Separation of Lifting-Line and Lifting Surface Velocities 

The difficulty of evaluating the integral of Eq.._ (29) (or the integrals 
of 04) (33) and ' (36) to infinity may be avoided by separating,, the trailing vor 
tex strength £Cr,«>into lifting- line and lifting-surface parts (Refs. 

9, 10, and 11). This is done by adding to and subtracting from Eq. (29) 
the quantity 

_l f' f $T jnm y e s * 5 ih 

z J 1 d ? W P 3 


The result may be expressed as 

, a r ( ?> R ^ • 

V (t \r.9)=\ j t 
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& 
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- It can be easily shown that vector V*(r,8) is equal- to 1 the - induced 
velocity V*(r) at radius r if the hydrodynamic advance coefficient Ai (f ) 

Xi 

is constant. For the case that is not constant, the approximation 
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will be made that V*{r,,e)can be replaced by V£(r) and considered to be 
constant as the point P moves along a helix. With this approximation, 
Eq. .(48) may be rewritten as 


Cm(Yj9) .. - _ U n M) 

Jr 2 + 4(r) ~ J l t*(r) 

where 
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Substituting Eqs. (43) , (45), and (74) into Eq. (73), we have 

&rVr) 

c m(Ke) " y*^) j \ir(UL%(r,<t>) + CL^dr,*) + 


-X x (r)(u^cr,<t>) -h U? ( * ? Ct ^ + u t 


(77) 


The integration with respect to for evaluation of the integrand of Eq. 
(77) is now limited within the bounds of ^(r) and # T (r). 

E-2. Acoustic Formulation 

The starting point of the acoustic formulation for propeller noise 
is an exact Ffowcs Williams-Hawkings governing equation of motion, which 
is a formal statement of the generalized basic equations of motion, con- 
cerning the sound generated by turbulence and by a surface in arbitrary 
motion (Refs. 16 and 23). This integral representation is referred 

to as the FW-H equation (18,1974). The FW-H equation shows that for a moving 
rigid body the acoustic density perturbation (f-j>) at a point X in the 
space-fixed coordinate system (x-system) and time t is given by 
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L R t- -|-*M 
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'T-r e 



aJ(z) 


di)(Z) 
TT = T e 


( 78 ) 


where 

a Q = speed of sound in undisturbed medium 
^ = instantaneous density 

1 = vector position of field point in the x-system 

Y = vector position of acoustic source point in the y-system 

t = initial time 
o 

t = non-dimensional observer time 
. R = vector distance between observer and source point 
R • = magnitude of R 
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= whole space at initial time t , 

S(t Q ) = moving surface 

0 C (t ). = volume enclosed hy S(t Q ) 

2 

T\ j = non-dimensional Lighthill stress tensor referred to 

V i = the ith velocity component of the moving acoustic source 

in the y- system 

V - (V.) - V +R-O.X z/v (for moving rigid body) 

1 o p p 

IX = angular velocity of z-system 

Z = vector position of acoustic source in the z-system 

V q = translational velocity of z-system 

M = V V/a 

p o 

a^ = the jth component of non-dimensional acceleration of the 

moving acoustic source in the y-system, non-dimensionalized 
with respect to V 
TT = dummy time variable 

T e = non-dimensional source time (retarded time) 

Mp = reference Mach number (tip Mach number of propeller) , 

- the jth component of force vector exerted on fluid by 
moving surface 

The first term of Eq. (78) shows that each moving element dl) (Z) out- 
side S(Z) is equivalent to a moving quadrupole source of strength dl) (Z) 
The second term shows that each element of surface area dS(Z) is equivalent 
to a moving dipole of strength -f^dA(Z) . The last two terms show that each 
moving volume element dl) (Z) , within S acts as if it emitted elementary 
waves which are the same as those emitted by a dipole source of strength 
-a^/Apd]) (Z) and a quadrupole of strength Vi Vj dl)(Z)and represent the 



42 


sound generated by the volume displacement effects. 

The consequence of the thin blade assumption is that the quadrupole 
sources T„ in Eq. (78) axe negligible and will be set to zero. The sur- 
face integrals may be replaced by a single-sided integral over the refe- 
rence surface, and the new source strength is the sum of the corresponding 
upper surface and lower surface sources. Furthermore, the integrands in 
the last two terms are evaluated at the reference surface and assumed to 
be independent of the n-coordinate, because the upper and lower surfaces 
are close together. 

Introducing the r-system into Eq. (78) and separating the loading 
and thickness effects, we have 


Acoustic pressure due to loading (force noise) 
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b. Acoustic pressure due to thickness (thickness noise) 
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c. Total acoustic pressure 
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It is noticed that the acoustic density perturbation has been replaced 
by the acoustic pressure, which is defined as the variation from atmos- 
pheric pressure. It is also noted that only steady sources are considered. 
Steady sources are those whose strength does not vary with time when 
viewed in the propeller-fixed coordinate system. 

A few remarks concerning various aspects of Eqs. (79) and (80) are 
in order. The notation £ S USed throu 8 hout this stud 7 to indi- 

cate that the quantity enclosed within the brackets is to be evaluated at 
source point ( <£ ) and retarded time X — f } *P) which is obtained 

by solving 

9 (a a x.f '<H = r e - 1 +M f |x- 7c*-., ?,•«>) | 

It is evident that Eq. (82) is an implicit equation for the required 
value of T e . If more than one solution to this equation exists (as it 
does at supersonic speed), each term in Eqs. (79) and (80) should be 
interpreted as a sum over all such solutions. At subsonic speed it has 
only one solution; for each source there is only one time at which it can 
transmit a signal to arrive at a given observer time, t. The next remark 

4 . ¥ — 

concerns the Doppler factor C * 1 - — • M, which occurs in the denominator 
of each term in Eqs . (79) and (80) . For supersonically moving sources, 
this factor can vanish at some point on the body and introduce a singu- 
larity with the resultant emission of Mach waves. 

The difficulty in evaluation of the integrals for arbitrary moving 
sources is associated with the determination of retarded time for a given 
observer time, t. This difficulty may be avoided by Fourier analysis as 
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is done in most work concerning rotational noise. Since we are concerned 
with the total acoustic pressure, Fourier analysis will not be applied. 
Instead, numerical iteration is required. 

Alternate forms of Eqs. (79} and (80) 

For the aeroacoustic study of propeller noise, the following alter- 
nates to Eqs. (79) and (80) are presented: 

Applying the chain rule to Eq. (82) shows that 
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Introducing the Doppler factor C , we have 
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Upon using Eqs. (82) and (84) to eliminate G from Eq. (83), we find 
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Hence, applying the chain rule to an arbitrary function f (X, T*e) 
shows that 
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derivative with respect to of a term of the form 
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where dependence on ^ and has been suppressed in order to simplify 
the notation. 

Upon tedious mathematical manipulation, it is found that 
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The terms inside the first bracket of each of Eqs. (87) and (88) 
represent the far field sound radiation from a point acoustic source while 
the rest of the terms represent the near field. It is interesting to 
note that the terms inside the third bracket of Eq. (88) dominate in the 
near field. 

Thus, taking the T derivatives under the integral signs of Eqs. (79) 
and (80) , we have the alternate forms 
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where 
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It is noticed that the dependence of acoustic pressure on , 

which -is the angular displacement of the y-axis with respect to the 
axis at t = 0, has been expressed explicitly in Eqs. (89), (90), and (91). 
Thus, for a given location, X, of the observer at time t, the instan- 
taneous acoustic pressure depends on the initial angular displacement Q q , 
of the propeller-fixed coordinate system. 
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E-3 Minimum Noise Criteria 

Since the initial angular displacement 8 0 of the propeller-fixed 
coordinate system relative to the space-fixed coordinate system cannot 
be controlled by the experimenter, it is treated as a random variable. 

For equally spaced and identical blades, the probability density function 
of this random variable may be expressed as 

f (0.) = a. 4 e. < a + air 

= q otherwise (94) 

where a is an arbitrary constant. In this case, the acoustic pressure 
(P “ P Q ) (X, ~t,0 o ) for a given X represents the entire family or ensemble 
of possible time histories which might have been the outcome of the same 
experiment. Since (p - p ) (X , ■£, Q 0 ) is a continuous function of the 
random variable d 0 > then (p - p ) [X, t, ) is a Borel function of d Q • 
From probability theory, the expectation of the random variable 
(p - p ) (X,- t, 6 0 ) is equal to the expectation of the function 

o — , 

(P - p ) (X, t, 0 O ) with respect to the random variable 6 0 

+ co 

E = j (f -%)(*,*,&) 16„ 

-CO 

which defines the mean square of (p - p Q ) (X, t, Q 0 ) . 

Substituting Eq. (94) into Eq. (95), we have 

Z'K -t A- 

Hv-r.fdw.) = ^j (t-r.fib-t, e.) d$ 0 
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Letting a = § - t and $ 0 = <p 0 ^ ^j+tand substituting into Eq. (96), 
we have 


2.7C 

t{[r-n){x,*,t>o)) =tj j(-p-n) (x>-t d to m 

0 


From Eq. (91) and all related definitions, it can easily be shown 
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It is easy to see from Eq. (98) that the mean square of the acoustic 
pressure for a stationary propeller is equal to its temporal mean square. 

Before proceeding further, it may be stated that Eqs. (25), (26), 

(27), (48) (or 77), (59), (65) (or 68), (6$(or 69), and (98) are the basic 
equations of the aerodynamics and acoustic (aeroacoustic) propeller problem. 
Most propeller problems are generated from this set of equations. 

A rather general problem of the design of a propeller for minimum 
noise is of finding the loading distribution and thickness distribution 
of a propeller blade that minimizes the ensemble mean square of the acous- 
tic pressure subject to constraints on the aerodynamic performance. From 
the calculus.- of variations, such a problem may be mathematically repre- 
sented by a nonlinear singular integro-differential Euler equation. So 
far, however, no one has succeeded in analyzing this problem. We shall 
consider the case for which the blade loading is modelled with a lifting 
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line and for which the thickness effects are ignored. 

Following are the equations for this simplified propeller model: 

Ideal thrust and power coefficients 
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Ensemble mean square of acoustic pressure (at t - 0, <£, = 0) 


_£rr 

E((-p-n)(.x,toM) = j Jp a*, m 

0 


It should be noticed that no loss in generality results from computing the 
ensemble mean square of acoustic pressure at t = 0 and 0 since the 
ensemble mean square of the acoustic pressure is independent of <§_ and 
the dependence of the acoustic pressure on time, t, may be replaced by 
the dependence on X and ($) . 
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PART 2. NUMERICAL FORMULATION OF 
OPTIMUM NOISE PROPELLER PROBLEM FOR THE 
SIMPLIFIED MODEL 


This part is concerned with the numerical formulation of the optimum 
noise propeller problem for the simplified model, and the computation of 
the lifting- line induced velocities. The study of aerodynamic problems of 
the propeller and of the "complex method" for constrained optimization had 
led to a nonlinear programming treatment of 'this model . We begin by expend- 
ing the aerodynamic and acoustic parameters in terms of Chebyshev polynomials 
and bivariate Chebyshev polynomials (31, 1973). 

The following new variables are introduced 


f(r) = 2r " ^ — 
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We note that 

p) =!,(') = i ; --i 

and 


S- r = 


(114) 


Since r(r) is continuous between the hub and the tip and vanishes at 


both ends, it may be approximated by a truncated expansion in Chebyshev 
polynomials 
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Let both the axial induction factor and the tangential induction factor be 
a PP ro;!C i raa Led by finite double Chebyshev series of degree N in both q and q Q 
of the forms 
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The interpolation is made over the points (q^ 
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are roots of T n+1 (q) = 0 and T n+1 (q Q ) = 0, respectively. That is 


cn+o 




_ow') r 12S + U 1 

1 TT W TD ) s = o 0) n 


(2s + or 


The double primes in. (117) and (118) indicate that the first term is h^ Q /4 and 

"t 2- "t cL "t SL * t 

h 00 /4 and that Ir ^ and h^.. and h^^ are to be taken as h^^/2 and h^^/2 

*£ 

and h^/2 for i > 0, j > 0, respectively. 

£ £ . . . V . 

The coefficients h^ and h^_. are determined by a biorthogonality proper- 
ty (31, 1973). 


7i n 

i i 

r=o s=o 


(n+0 Cn+O 


Cvw) c*+0 


2 2 Tlj (%y 5 ^ ( tr i ?os ) 


= ( Y\-t Q 
4 


= (n-H) 

2 

I 

— ('TH-l) 

- 0 


or 


X = K * 0 

J » * * 0 

4 = R 4(3} j - 1-0 

A. - K =0 , j =r J> =40 

i « K , A 
■i 4 K J -i. 


so that 


' ft 


4- 71 u 

CH+0*£o|* 


2* 

^0 


(n+Q (nti) C***') (n+» 

(, ?r ’ ^ V ? tos ^ 


( 119 ) 



58 



A backwards recurrence formula for the evaluation of the sums of Eqs 
(119) and (120) is given in Appendix B. 


A. The Evaluation of Induced Velocities 

Substituting Eqs. (117), (113), (113) and (116) into Eqs. (SO) and (51) 
respectively, we have 


hp . IWTWv 1 


(121) 




+i 


TwC^c) J9 (122) 

FTft-p 


Furthermore, applying the solution for a Cauchy singular integral of the 
form 


l (VPfFx 


= V l)*M (?) 


we obtain 


ut (?) = I { 1 1” TK?)Ti tp U*-, «) 

H irt=i 1 % V t=o^-0 d U 

+ 11 ” hi TAWy (!•_,(?) I 

X=o^m-H 0 d • J 


2N+MH 

= Z * Cf) 

4.-0 


( 123 ) 
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<‘i>- 5, i £ f: T ‘‘i> !<!' u«<() 

+ ££” &TK»"L(fl'U«(f)| 

£=oj = w+i <j Q * 


2N+M-1 


= £ ' U tji T^p 

* —rt * 


(124) 


where u . and or . are the coefficients of the Chebyshev expansions of the 
a,i t,i • • 

axial and tangential components of the induced velocity , respectively. 


B. The Evaluation of Hydrodynamic Advance -Coefficients 


Combining Eqs (123) , (124) and (113) , we have 

2N+H-1 

w? \= «»>(%■» 5.’ 

A 0 2*/+M-i 


(125) 


nv 

A-p 


’ T Vt :Ti(^ 


It should be noticed that Eq. (125) is an implicit equation for the required 
value of X^ (q) for a given T (q) since the induction factors are functions 
of A^(q). Equation (125) will be used to obtain the new approximation for 
X i (q) from the current value 1 X^(q) in the iterative scheme for solving the 
non-optimum propeller problem (8, 1952) to obtain an initial solution for 
the input of the. complex method. 


C. The Chebyshev Coefficients- of Circulation 

Let q^ ^ be the roots of 1^ (q) = 0 and X^ ^ the values of X^(q) at 
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Using the Eq. (101) , and Equations. (.123) and'.C124)., we obtain a 


system o£ linear equations for M unknowns, G^. 


CM) 

At 

Xp 


M r N m o, , f ' L u . . (M) 

- * '£,(£?< 1 V 


CM) - -t 
2 Aji hli 


N M„ «. 

+ ir (h* 


CM) * 

1 ni /U\ V 




& 


J? = I CO M (126) 


3. "t * 

where lu ^ and lu.. are obtained from Eqs.. (119) and (120), respectively,. 


and 


CM) 


(M) 




(M) 


= WTjCt, ) U-,(|, ) 


(127) 


The coefficients are determined by solving Eqs. (126). 

D. The-. Evaluation of Thrust^ and Power Coefficients 

Using Eqs (99) and (100) and the orthogonality property of the Chebyshev 
Polynomials, we obtain expressions for the ideal thrust and power coefficients 



B(i-ft.)7r I it. A. q 
z l 


H- 


i r ft. 

2- A-p 


M 

+ XI Gw ( 

m=i 


U 


t,m+ 1 



G 


a 


(128) 
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c p^= - b t x jr x W + G: 


0 -*A) 


Vf 


w*=l 2! ^ r ' irl ( 


M 


+ 


‘‘fc 


Tn=| 


jw-21 " Ua, w+ij 


(129) 


E. The Evaluation of Acoustic .Pressure 

To compute the instantaneous acoustic pressure, we separate the kernel 
K(X, (JJ) , p, <j> Q ) of-Eq. £102). into three -parts 


= K, (X,®, •?,<?<,) 


t )<£(*) t K t (X,®, f,vu*(f) 


(130) 


where 


ko(x,®, ?,&)= E 

K-l 


_P 

R 


? ' c « {n f Vf ( 

|->*®c«(4> l+ £ K+Mf . i ^ + J_ 

j- (-£• ® - h f ) + + + 


t M, 


)( 


( 131 ) 
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B *P r 

4. J_ ^_!. J1_ + ) -5-y^ Mfp)|l (132) 

' J T=r e 


V*.»,!,t) = g (trc“ 




i-m‘ 


4- 




(133) 


^r=r e 


Furthermore, let K o ,(X,$, q Q , <j> Q ), KJX,®, q Q , tf> a ) , and KJX,®, q Q , 4> q ) 
be approximated by finite Chebyshev series of degree L in q Q of the form 

l , 

k„ U.®.f.»*.) = Z P.,(cCJf,M.)Ti{(?.) a 34 ) 


*=o 


L(XA%,P.) = Tk«o) 


*=0 

L 




(135) 


(136) 


where P Qjk (X,<g>, <f> Q ) , p ajk CX,<3>, * 0 ), and P k (X, ©,' <J> Q ) 
can easily be obtained by using Eq. (B-2) of Appendix B. 

The advantage of separation of the kernel is that the coefficients P 

O jK 

(X,®, <f> 0 ) > p a k (X,(g), d> o ) , and k (X,(g), $ o ) may be computed once for all 
for given B, M , X,(g), and <#> o since they are independent of the helixal 

vortex system behind the propeller. 

Substituting Eqs. (130) and (115) together with Eqs. (134), (135), 
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and (136) into Eq. (102) we. have 


M L 

Q'm 

K s 


2 4 


m=t 


ML, . 

B,m+0 + 21 -Gm 21 (fajK^m^K + ft,K Et,m,i< )f 
m=i k=o J 


where 


+( 


P*,« )K = j fFjJ" T K a 0 ) llm-cp uUt) d|, 

-1 


2L 

a 


2I/+MH 

21 ( fy° ^ K>0 ^ Sjj K+m-l 

J — 0 


+ £j,0 £&c,Tn-i + jx-m-Hj - w + K+l 


ft, 


w»,K — 


~ Sjj,0 “ Sjy (.K-W -t |) 

= j j i ~ TxC^j,) (Jw-i ($ 0 ^ MV> 


2N+M-1 


g 21 3 U tji ^ %o ^K,o £*1-1 + ^',K+ W '“ 

j s ° * 

+ £j,o S K/'m-t ■+ <£^‘, J jc-m-ti 1 “ » vn+ K + 1 

- £j,o 


$A,j =• O if i r j 

= | if i = j 


(137) 


(138) 


( 139 ) 
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F. The Evaluation of the Ensemble Mean Square of Acoustic Pressure 

We start by writing Eq. (112) as ^ 

0 

H 

~~Z J(^~ {X > ® j *£) d’fj (140) 

-i 

where 

+. = \ ( i + 1 ) 

Now, let n • and W. be the abscissas and weights of- the K -point Gauss-Legendre 

•) J 

formula. Then Eq. (140) may be approximated by 

K 

E(£f-?n(X>®> +•>) = |lWj C T-rJ(X>®> ?j) (141) 

where (p - P Q )(X,$), Tj..)' are evaluated using Eq. (137). 

It should be noticed that the retarded time and the distance R are com- 
puted by using Newton's method (see Appendix C) . 

G. The Nonlinear Programming for the Simplified Propeller Model 

In this section we are concerned with the formulation of a nonlinear 
programming model for the simplified propeller. To begin, we assume that 
the number of propeller blades, B, the advance (or forward) Mach number, M_, 

r 

the tip Mach number, M , the distance between the observer and the center of 
the propeller, X, and the azimuth angle, $), of the observer are known. 

Investigating the evaluation of the aerodynamic and acoustic quantities 
of the propeller, we find that for a given A^q^, a H these quantities are 
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determined. That is, for a given configuration of the helical vortex system 
behind the propeller all aerodynamic and acoustic characteristics are fixed. 
This suggests that it is possible to find a configuration which satisfies all 
the specified constraints and produces minimum noise. A nonlinear programm- 
ing model is established to facilitate numerical determination of this opti- 
mum configuration. 

Let he the values of A. (q ) at the zeros of T t , (q ) . A nonlin- 

j l hr J+l vn o 

ear programming model for the simplified propeller has the following form: 


Maximize 


e*(™, 


< 7+0 

A z , 




Subject to 


crw) 


- A Ul j } = 1 0) Ttl 


-TL 


4 c Ti z c. 


TU 


'2L 


/ C B i 4 C 2 U 


(142) 


where the ideal thrust and power coefficients are. regarded as implicit vari- 
ables while ^(J + l) ^C J+ i) are t * ie ex Plicit independent variables. 

1 ’2 * J+l 

The upper and lower constraints are either constants or functions of the in- 
dependent variables. It is noticed that the value of X^(q Q ) at any point q Q 
is interpolated from the polynomial interpolation of degree J which exactly 
fits A i (q Q ) at q/ J+1 ' ) , j = 1, 2, ... , J+l. 

The algorithm of J. A. Richardson and J. L. Kuester (32, 1973) based 


on the "complex" method of M. J. Box (30, 1965) has been modified with two 
feasible starting points as input to solve the nonlinear programming (142) . 

The constrained complex method is a sequential search technique.- Since the initial 
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set of points is randomly scattered throughout the feasible region, the pro- 
cedure should tend to find the global maximum. The first feasible starting 
point is the solution of the aerodynamic optimum propeller. The second fea- 
sible starting point is any solution for a non-optimum propeller. These two 
feasible solutions are generated by using the iterative technique given in 
Ref. 8. The procedures are described below. 

A propeller having a constant hydrodynamic, coefficient is called an 
aerodynamic optimum propeller since its ideal efficiency is, according to Betz, 
the greatest that can be obtained for a given propeller advance coefficient 
and ideal thrust coefficient. 

Procedure for Finding Aerodynamic Optimum Solution 

1. Specify the loading coefficient (or C p ^) which satisfies 
the implicit constraint of Eq. (142) 

2. Assume A 0 ^ = A. , % = £(1)M 

■ 3. Solve the system of Eqs. (126) for 
4. Compute (or C p ^) from Eq. (128) (or Eq. 129). * 

Repeating steps (2) through (4) for several values of A^, the- dependence 

•k 

of the loading coefficients on A^ is derived from which the proper value A 
is interpolated. Thus, the first feasible solution is obtained by setting 

Aj = A. 2 = 10) T"ti (143) 

Procedure for Finding a Non-Optimum Solution 

1. Specify the loading coefficient C Ti (or C pi ) which 
satifies the implicit constraint of Eq. (142) 
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2. Specify a characterising function for the circulation 

3 = z ^ u« ({) (1443 

3. Relate the circulation r(q) which satisfies the loading coefficient 
to the given characterising function by a factor k which is indepen- 
dent of q 

M 

rep = m = ?«iL-icp (ms) 

“ m-l 

or 

Gw ~ ^ 

4. Assume = A p j = 1(1) J+ 1 

5. Compute u . and u . from Eqs.. (123) and (124), respectively, by 

3.J 1 "C y 1 

replacing G ffi by g ffi 


lu cp = 

2N+M-1 

Wo«.,xTiC%) 

-t= o 

(146) 

= 

fc 21’ 

(147) 


A.-0 


6. Substitute Eqs. (145), (146), and (147) into Eq. (128) (or 129) 



C-u = R C T , + 

K Ct2 

(148) 

or 

Qi = * c n + 

K*C r t 

(149) 


where 
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cu (},<!+»*) 


z x 


"T 2. 


S - * ~ '»»+*) 




c 2{ = **i T " U< t4 ^> 4 8* 


Cp 2 — r^~ | ( l * t'k) %• ( W-A.jW-1 *“ W-h) 

l* 4* v m=l 

M 

E 

7»=t 

7. Solve Eq. (148) (ox 149) for k. 

8.. Obtain a new approximation set of xj J+1 ^ by substituting Eqs. (146) 
and (147) into Eq. (101) 

Repeat steps (5) through (8) until a satisfactory feasible solution is 
obtained. 

The numerical computation* has been programmed for the CDC Cyber 175 com- 
puter. The FORTRAN listing for the program may be found in Appendix D. De- 
tails for use of the program is outlined in the main program. In addition to 
the nonlinear program Eq. (142), this program is capable of solving the follow' 
ing three problems: 





Problem X: 

Given: B, r^, V p , X^ (or Mp and M^) , X i (r) = X^> X, and 

Determine: Aerodynamic optimum circulation distribution. 

Induced velocity components 

Si* C Pi 
Ideal efficiency 

Ensemble mean square of the acoustic pressure (optional) . 
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Problem II: 

Given: B, V p , X p (or M p and M p ) , X^r) = X ± , 

C Ti (or C pi ), X, and (g) 

Determine: Aerodynamic optimum circulation distribution. 

Induced velocity components. 

Ideal efficiency ' - 

Ensemble mean square of the acoustic 
pressure (optional) 

Problem III: 

Given: B, r^, V p , X^ (or Mp and M^), 

C Ti (or C p .), X, and (j?) 

The type of characterising function of circulation 
Determine: Circulation distribution. 

Induced velocity components. 

Ideal efficiency 

■ Ensemble mean square of the acoustic 
pressure (optional)) 


Results of sample, calculations are present in Part 3 . 

Having determined the acoustic optimum circulation and the hydrodynamic 
.dvance coefficient distributions for the propeller, the actual shape of the 
•lades remains to be determined by lifting- surface technique. It is only 
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necessary to select an appropriate chordwise circulation distribution and a 
thickness form at each blade section to evaluate the lifting-surface veloci- 
ties. The mean line at each blade section is then obtained from Eq. (48). 

The reader is referred to Refs. 10, and 35..for detailed numerical - computation . 
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PART 3. APPLICATIONS AND NUMERICAL RESULTS 

The main purpose of the computing program based on the technique develop- 
ed in Part 2 is to solve the nonlinear programming Eq. (142) to find the acous- 
tic optimum circulation within prescribed aerodynamic constraints- In order 
to give some verification of the present technique, some computed examples are 
given in the following 

Numerical Examples 

1. Given: B =5, r, = 0.2, A = 0.19966, V c = 1, and A. = 0.27211. 

1 h 5 p F l 

Determine: Aerodynamic optimum circulation distribution, 

induced velocity 
thrust and power coefficients, 
ideal efficiency 

The results for M = 10 and N = 10 are shown in Table 1. 

The last three columns in Table 1 show a comparison of results for the 
same propeller obtained from Ref. 10. The agreement between these methods 
supports the validity of the present aerodynamic model and computing program. 
Finally, the thrust and power coefficients have been determined, = 1.23324, 
Cp^ = 1.68071. From these, the ideal efficiency is 0.73376. From the well- 
known relation for an optimum propeller, ri^ = A^ Vp/A^, there is obtained 
n. = 0.73375. 

l 
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Table 1 

Results for an Optimum Free-Running 
Five-Bladed Propeller with. A = 0.19966, Vp = 1, A.= 0.27211 



G 1 ~ 

0.03301036 

G 2 “ 

0.0034127-2 

G„ * 
o 

0.00125070 


G 4 = 

0.00018100 

G S = 

-0.00005210 

G 6 « 

-0.00001478 


G 7 = 

-0.00000722 

G 8 = 

0.00000089 

G 9 = 

0.00000028- 


G 10 = 

0.00000027 






C Ti = 

1 j 23324 






C Pi = 

1.68071 






\ = 

0.73376 


- 

Kerwin f s 

Lerb * s 
induction 





vortex- line 

factor 




* 


method 

method 

r 

q 

u* 

a 

U t 

F • 

F 

r 

0.2 

-1.00 

0.12674 

-0.17197 

0.0 • 

0.0 

0.0 

0.3 

-0.75 

O’. 19899 

-0.18047 

0.01945 

0.0197 

0.0196 

0.4 

-0.05 

0.24802 

-0.16869 

0.02582 

0.0260 

0.0258 

0.5 

-0.25 

0.27992 

-0.15233 

0.02955 

0.0297 

0.0296 

0.6 

. 0.00 

0.30098 

-0.13650 

0.03171 

0.0318 

0.0317 

0.7 

0.25 

0.31522 

-0.12253 

0.03252 

0.0325 

0.0325 

0.8 

0.50 

0.32521 

-0.11062 

0.03142 

0.0314 

0.0314 

0.9 

0.75 

0.33240 

-0.10050 

0.02634 

0.0262 

0.0263 

1.0 

1.00 

0.33719 

-0.09154 

0.0 

0.0 

0.0 
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2. Given: B = 2, r, = 0.2, V„ = 1, A = 0.26932, M„ = 0.2, 

h F p F 

M p = 0.7426, X = 4.272, and ( 3 ) « 1,212' rad. 

Constraints 

o.ZS i Jlj ^ 0.5 

0.74174 C-ri 

Determine: Acoustic optimum circulation distribution, 

induced velocity 
thrust and power coefficients 
ideal efficiency 

For this particular case only the feasible aerodynamic optimum solution 
is necessary to be used as the initial- solution in the complex method. All 
points, which were generated from random numbers and constraints, converged 
to A^ » 0,28 quite rapidly. The results are shown in Figs, •-[9’) through (11). 
It is not surprising that the acoustic optimum solution is the aerodynamic 
optimum solution that has minimum thrust coefficient since all loadings in 
this case are similar as shown in Fig. 10. This provides a simple numerical 
check of the technique developed in Part 2. 

It is noted that Fig. 11 is plotted based on the following definition of 
the total sound pressure level: ~ 

Et(p-p 0 n 

Total Sound Pressure Level - 10 log 1Q 9 

p ref 

2 p o V r» 2 

= 10 log 1Q E((p-p o r) + 20 log 10 dB 

In order to apply the present technique to a more general case, the 
third example is chosen as follows 
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0 0.2 . 0.4 0.6 0.8 1.0 

r 

Fig. 10 Relation between aerodynamic optimum 


circulation and X- 
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3. Given: B = 2, R = 1.524 m = 0.2, V'= 67,056 m/s 

P P 

» = 163.363 rad/s, V = 1, 

F 

Total Thrust = 7116.8--N 


Constraints 

0 . 15 xj J+1) <_ 0.6 j » 1(1)10 

0.35560 < C Ti 

Determine: Acoustic optimum circulation distribution, 

induced velocity, 
thrust and power coefficients 
ideal efficiency. 

In calculation of the ensemble mean square of the acoustic pressure, 
the following values were used 

, , 3 

air density p Q = 1.22 kg/m 


speed of sound a Q = 335.28 m/s 


In the. numerical computation, J, L, M, and N were taken to be 9, 10, 10, 
and 10 respectively. The two feasible solutions corresponding to = 0.3556 
are shown in Fig. 12. The characterising function is also shown in Fig. 12 ’ 


g(q) = 0.03267 - -0.00083 tt^q) -0.00467 U 2 (q) 
- 0.00243 U 3 (q) -0.00027 U 4 (7) 


After 19 iterations the factor k was found to-be 1.01916 and all X. 

3 

satisfied the convergence criterion 


CJ+D 



0 -. 0-6 



B = 2 , 

r h“ °' 2, 

< 

~n 

li 

p— 0-2, 



M p = 0.7 4 2 6 , X = 

= 4.27 2, 



0.0 5 

__ 0 =1.212 rad 





C Tj = 0.3 5 5 6 

r acou. 

opt. 




k = 1.0191 6 

/ \ 



0.0 4 
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0.0 3 
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// 

warn 
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mH 






nil 






HH 

0.01 
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0 - 0.2 0.4 ^ ' 0.6 0.8 1.0 


Fig. 1 2 Acoustic optimum and feasible 


circulations. 
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x? J+1 \ - xP +13 < 10 - 6 

3 , n+1 j,n - 


The adoustic optimum solution, was ' found ;after.-117' iterations , Convergence - 

for the= complex method was assumed when the objective function values at each 

-8 

point were within 10 units for 4 consecutive iterations; The toral execu- 
tion time was 164.531 central processor, seconds. 

The ensemble mean square pressures for the acoustic optimum solution, the 
feasible aerodynamic optimum and the non-optimum solution were found to be 
0.1082xl0~ 4 , 0. 1166x1 0 -4 ,. and 0.1132xl0~ 4 , respectively. The corresponding 

2 v 2 

root mean square of the acoustic pressures were 18,044 N/m , 18,752 N/m , 

2 

and 1*8. 456 .. N/m „ The computed results are' shown in Table 2. 

To give a comparison of the order of magnitudes of the root mean square 
of acoustic pressures, we note that the root mean square of the pressure for 
the fundamental harmonic of the propeller with the same operating conditions, 
except slight difference in power coefficient is 15;-4 N/m obtained from 
Fig. 5 of Ref. 13. It should be noted that no far field or near field assump- 
tion is made in the present formulation while the root mean square pressures 
shown in Fig. 5 of Ref. 13 were obtained with the assumptions that the field 
point is in the far field and the radial integrals are replaced by an effec- 
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Table 2 


Acoustic Optimum Solution for a Free-Running 
Two-Bladed Propeller with Operating Conditions Given in Example 3 


G 1 ■ 

0.03651389 

' G 2 

= -0.00863877 

G 3 

= -0.00565379 

II 

U3 

0.00087471 

G 5 

= 0.00018965 

G 6 

= -0.00171678 

II 

f- 

-0.00000712 

G 8 

= 0.00059500. 

G 9 

= -0.00010552 

G in .= 

0.00031216 






C„. = 0.35560 

Ti 

C^.. = 0.44604 

Pi 


.79724 


Ensemble mean 

square = 0, 

1082X10" 4 

r 

q 

u* 

a 

U t 

r 

0.2 

-1.00 

0.08376’ 

-0.00831 

0.0 

0.3 

-0.75 

0.09020 

-0.08432 

0.02606 

0.4 

-0.50 

0.21581 

-0.15674 

0.03961 

0.5 

-0.25 

0.28651 

-0.18653 

0.04639 

0.6 

0.00 

0.23990 

-0.12160 

• 0.04225 

0.7 

0.25 

0.15667 

-0.06066 

0.03277 

0.8 

0.50 

0.11000 

-0.03649 

0.02328 

0.9 

0.75 

0.07638 

-0.01579 

0.01268 

1.0 

1.00 

-0.12592 

0.03202 

0.0 
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APPENDIX A 

EVALUATION OF INDUCTION FACTORS 

Since we axe only interested in the axial induction factor and the 
tangential induction factor, the radial induction factor will not be considered. 

Integrating by parts, the tangential induction factor, Eq. (55), may be 
expressed in terms of the axial induction factor as 

With this relationship, only the evaluation of the axial induction fac- 
tor needs to be discussed. In. addition to Wrench's modified formulas (28, 

1957), an alternate method is presented. 

\ 

I. Wrench's Modified Formulas 

The Wrench's modified formulas for evaluation of the axial and tangen- 
tial induction factors may be summarized as 


1*0;. f) = e (i-^H i -2B3.F,) C r < f ) 


= 26 2 «(i- -y-)F 4 


( r>f ) (A-D 
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lt(r,tt=- 2 BH.Ci- r^JF, ( r < J ) 

= -B(l-^-)(l+2B1f.Fa) (OJ) 

where , 

^ = *IcF) 

H- =, -L- 

/o Ai(f) 

c ] (HjLf {_[ + i f ??*+ 2 

r ‘~ 2B«fc V 1+ / | wT'-j + 246 V (i + tffz 


+ 


5 2 3 -2 

( 1 + ?*)** 



( + - 


J 

JT’-I 



r ^ i i i ± it \ » r i » / f^-^2 

**• 26^1 i+^ z ; l U-l 243 V {{+ ^Va 


+ 




/ JiTlT - 1 P fv l+ ~ J /- <^ s 

v ? JTTtf-iJ e 


(A- 2) 


CA-3) 


CA-4) 


(A-5) 
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Fox detailed derivation of Eqs. (A-l) and (A~2')i the reader is referred 
to Ref. 28. 


II. Alternate Method 


For numerical computation, the integral of Eq. (54) is split up into 
four parts: 


I a (r,j) = + Ia(r.?) + I*fr.?) + ifo?) (A ' 6) 


where 


2?r L 

T'(hj) = f (f- r)Uf-r<*«A) ; du. 

* i { f ’ +t >- 2S ^ 


(A-7) 


27ri- 

T a o^ = ( f <dzl2Jll + dJu . 

* ' J *= 2 [f+r'-sr? Sc) " 


e 


T 3 (re) = f c s-r) $ C ? - r 
ft J { r^-zf r cAACd-Ajcf )^ 1 ] 3 / 2 



(A-8) 


(A-9) 


T 4 (re)= f y ~- ( JzlllAl L Sk)) dM. 

* ' tlF' f?*+ r x -2?rac^+WHr^cs)X} % ^ 


(A- 10) 


where L is an integer, chosen such that 
2v L X i (p) > 2 

j 

and 0 < e < < 1 
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1. Evaluation of I*(r,p) and (r,p) 

3- 3. 

To evaluate- these integrals the interval of each integration is divided 
into several subintervals. A 25-point Legendre-Gauss formula is 'used in each 
subinterval. 


2. Evaluation of I (r,p) 

3 


In this case, y is small, and cosy may be substituted by 1 
With this substitution Eq. (A- 9) may be approximated as 

,6 


u 2 /2. 




Cf-r) S {(f-r) + 

{ (S-rf + Crf + Alc?)^ 1 ] 


3 /z 


CA-11) 


Upon integrating Eq. (A- 11}, we have 

Ia(r.f)= fe 


jsr(f-r) 

■J ( ?-rj 2 +(rJ +>!(?)) 6 ^ (Irf 


“t* 


I 


. ^ ^ f 

JCf-r^+Cfj+A^e 2 ' JT7+^5T?r o ^ 


6 i ry + fcy-r) 2 -^ (rg -hAffi) 6* 

|f- M 


) 


CA- 12) 


3. Evaluation of I (r,p) 

3 


Since r and p are no greater than one and 2ir L A^(p)>2, the integrand 
of Eq. (A-10) can be expanded as a hypergeometric series F Q (32, 1969} 
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= (?-r) Z A«Z C2t n 'V a { f*COSF(u,M+J,jL) 

11=0 4=0 


- ?rcosF(u,2n+3 > 0} 


(A-13) 


where 


p * 

q = 


A = 
n 


Cr 2 + p 2 )/4 


- rp/2 ’ 2 

re!) rc™-o (Mf)u)- 


u 


* L 


COSF (u,n,i) 



8 i 

22 C + Ut) 



dt 


CA-14) 


4 

Once COSF have been calculated, I a (r,p) can be evaluated by direct sub- 
stitution- It is noted that COSF does not depend on (p) but does depend on 
the number of propeller blades, B, and parameters u, n, and i. 
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Evaluation of COSF 


The following are some useful formulas for evaluating COSF. 


B 


T" 1 (fc-t) = 0 for any n and B 

t-l S 


B 277* 

y i2L(K-0 

K-U 5 


_a>c 




— B if n/B is an integer 

c 0 if n/B is not an integer 

* 

A, 


<2&<L $ =. — T- -4* 


C • | 

__ ±: 4 . — • y C/_4 

2 1 " 2 l>1 fe & 


5VC-1 
OSKt & 


1 X 3^-1 

73 * z c H ^^r oe 

2 3=1 3 - 


Using these formulas with the sine-cosine product relations, we have 

D „ 1 A* 'Xp. . — ~ 

Z Ut) =B(-p. + y £ Q-j ^ 2 5 ut ‘% 8 ) 

K=l 4 4 3“* 0 


B 


E ^ Wlff-Ut) = ~ (E c/j Ut 5 pg-o/g 


where 



0 

l 


if x is not an integer 
if x is an integer 


(A-15) 


(A- 16) 


(A- 17) 


(A-18) 


{A- 19} 


(A- 20) 
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Now, with. Eqs. (A-19) and (A-20), we are able to express COSF in terms 
of summations of the generalized cosine integral. Cl. Dividing Eqs. (A-19) 
and (A-20) by t n and then integrating from 1 to infinity with respect to t. 


we have 

C0SF(U i r\ 1 2l) 


co 2> 2 .X 

Z ( SicfUt ) 

— i d-t 

t * 




= B { + 7 X -1 |r, Cl >•») J 


(A-21) 


CQ ^5 3^-1 


COSF(U,T»,U-0 = 


± t+C C^K + U-fc) 




it 


t 


ru 


4rr £ C H- Cl ( l 1 }-') a »«) 

2 3=t <r 


(A- 2 2) 


where Cl is the generalized cosine integral defined as 


.GO 


CIO, 71 ) = 




cdz 


n. 


4± 


t 


Asymptotic expansion of Cl (a,n) 


(A- 23) 


Applying integration by parts to Eq. (A-23) , we obtain (33, 1957) 
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Cic a , 27 l) = £ (ft) + 4 57 <? 0 <>) +>(«)} (A- 24 


Ci(a j2W )=|^i^Q:c a .)-^^) + C*WJ 


(A- 25) 


where 


£0 


^.Z(a)=-j -i^d-t =-| + J.(a) 


A. 


CL 

Si C.0-) = j ‘‘‘t 


2? i f s&n'L ?4- 
• I + j ar (A- 26) 


(A-27) 


Ci(ft)= f i2±f. d-t 

J TL 

CO 




(A- 28) 


(A- 29) 


'(A-30) 


The asymptotic representations of Ci(a) and — - Si (a) may be found in 
those books which deal with special functions. The following are taken from 
Ref. 34, with some change in notation 


aw 


0 L, 






CA-31) 



Coo c< - 
<X* 


/?%> + 





(A-32) 
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where 



<:-»*(>*) ! 


-f Oc 


— 2rt“X 

|o.l 


) 



ei/fty+oj 


+ 0 ( 


-an -3 

1*1 ) 


Using these expansions in Eqs. (A- 24) and (A- 25) , we obtain 

y. JJT-l 

ctcv«) - ftca) + 


CICa^n -0 


rt 2 «~ z 

(rO f o^c^ n 

CMI-20! t ^ ^ 



where 


P,»> = Pn^Co) 

a B w = Q*o) 


For convenience of numerical computation, "we define 



J_f J_ q»)(an-H) f. 

M 1 a 2 - V 




a 


Czn+4-) frn+s,) 

a. 2 * + 



C*V 4 \ -i" / 1 _ &Si Otfiw) 

5 '" < - J ~aM <V 



<>w*3)£ngO 

aF 



£2Y\+5) 0*^+6) -5) J 


(A- 33) 


- CA-34) 


(A-35) 


(A-36) 


(A-37) 


(A-38) 
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tfc) (1- fy -»-)) ] (A-39) 


i ] ]_ 2nferw) /, foa*a)(*m-3) { | x\ 7 (A-40) 

AM A 2 * V 0.* ' A 1 v j 


Then, 

^ ^Sr 

Cl OMtO = (a) oi^o ~ !i m (a) _M*s o o (A _ 41} 

ci (A, 2n- i ) = []*(*) - T^>) ^ (A ' 42) 


With a proper choice of L and e, this alternate method will yield induc- 
tion factors to any desired accuracy. It is noted that this method can he 
extended to the evaluation of both the tangential and radial induction factors. 
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APPENDIX B 


EVALUATION OF COEFFICIENTS h*. AND h*. 

ij 13 


ci t 

To evaluate tlie coefficients and lu.., we give two backward recur- 
rence formulas for evaluating polynomials in Chebyshev form 

a ~r r *\ _ 




where 


— 3<^+2. ^ 


3-k — Bk+i - * +4^ K - i\, 01 - 1 , •** jO (b-i) 


2 * 

= % ( Be" 


where 


£, k = 2 ± Bs+\ - 8 K+2- + Cn=+> 

t =n i - 1 


To evaluate the coefficients h. . , we put 

ij * 

-r l4 tm0 

Cpg — la (. 15 *- * Qos ' 


£=m,m-V" » 0 


(B-2) 


CB-3) 


and 


, -Cn+O \ 

P^ = I e C rs T^(t s ) 


so that Eq. (119) becomes 

-n 


- A 


CB-4) . 


CB-5) 
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Define 



= 


2 ( 71 + 1 ) 


-Then, 


_r r q \ ~r ( *y ^ +,i \ 

TAlos ) ~ ‘iS-H C ) 


and 


tic c +0 


) — Tir+i ( ^ 


Cmi) 


) 


Upon substituting Eqs. (B-6) and (B-7) into Eqs. (B-4) and (B— 5) , 
spectively, we have 



~ [2 


2 Sir \ 




4 

(ntr) x 


n 


r=o 


Tzr-ti 



) 


The advantage of Eqs. (B-8) and (B-9) is that and can be 
evaluated by using the backward recurrence formula (B-2) 

The coefficients h?.. are obtained by replacing I a (q^ n ^ ” 

in Eq. (B-3) by I t Cq r Cn+1} , qj n+1) ). 


l os 


l os 


(B-6) 


(B-7) 


re- 


(B-8) 


(B-9) 


C n+ l) 



APPENDIX C 


EVALUATION OF RETARDED DISTANCE R 

= T- tV {- ® 

+ J Mp «-/<® -t- Cl - ') 0+-^"- 2 ( c -15 







(C-2) 


£o “ ‘ppjjp ® ■+ J l- MjAj*? © } 

Then, R is found by the Newton method: 


■f ( Pn> >) 


— ^n-i 


(C-4) 
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C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


APPENDIX D 

FORTRAN LISTING OF COMPUTES PROGRAM 
PROGRAM PROPEL (INPUT, OUTPUT, COSGRA , PK0ATD, 

$ TAPE1 =COSGRA, TAP22=PK0 ATD,T APE6=OUTPUT) 

PURPOSE 

TO SOLVE PROBLEMS I THROUGH IV LISTED BELOW 
BASED ON THE' TECHNIQUE DEVELOPED IN PART 2. 7 

I 

DEFINITIONS: 

VP=REFEBENCE VELOCITY {USUALLY, 

VP= ADVANCE SPEED OF • PROPELLER) 
NBLADE=NUMBER OF PROPELLER BLADES 

RP= REFERENCE LENGTH {RADIUS OF PROPELLER) 

T (N ,Q) =CHEBYSHEV POLYNOMIAL OF THE FIRST KIND 
OF DEGREE N (-1.GE, Q ,LE.1) 

U (N,Q)=CHEBYSHEV POLYNOMIAL OF THE SECOND KIND 
OF DEGREE N 

RH=HUB RADIUS OF PROPELLER 
RC=RADXAL COORDINATE OF BLADE SECTION 
RH ,gt,:rc ,LE. RP 
OMEGA= ANGULAR VELOCITY OF PROPELLER 
RAMDAI=HYDRODYNAMIC ADVANCE COEFFICIENT 

Q={2.*R-HUB-1) /(1-HUB) {-1,GE, Q «LE,1) 

WHERE HUB=RH/RP , HUB '« GE* 7 R=SC/RP . LE. . 1) 
T=TOTAL THRUST 
P= TOTAL POWER 

MSP=ENSEMBLE MEAN SQUARE OF ACOUSTIC PRESSURE 
RAMDAP= REFERENCE ADVANCE COEFFICIENT 
(RAMDAP=VP/ (RP*OMEGA) ) 

VF=ADVANC2 SPEED OF PSOPELLER/VP 
FM=ADVANCE MACH NUMBER OF PROPELLER 
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C 

c . 

c 

c 

c 

c 

c 

c 

c 

C' 

c 

c 

c . 

c 

c 

c 

c • 

c 

c 

c 

c 

C' 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TH=TIP MACH NUMBER 

CT=THRUST COEFFICIENT—!/ (Q» 5*DENSITY* 
VP**2*PAI*RP**2) 

CP=POWER COEFFICIENT=P/(0»5*DENSITY* 
yp**3*p AI * R p**2) 

ETAI=IDEAL EFFICIENCY=CT* VF/CP 
ENSMSQ= NON -DIMENS IONAL ENSEMBLE MEAN SQUARE 
OF THE ACOUSTIC PRESSURE 
=MSP/(DENSITY*VP**2/ (4*PAI) ) 

✓ 

XDISTA=NON-DIMENSIONAL DISTANCE BETWEEN THE 
OBSERVER AND THE CENTER OF PROPELLER 
AZIHUT=AZIMUTH ANGLE OF OBSERVER (IN DEGREE) 

(C .GE. AZIMUT • LE* 90 (DEGREE) WHEN 
THE OBSERVER IS BEHIND THE PROPELLER 
DISK 

UA=NON-DI MENS IONAL AXIAL COMPONENT OF 

INDUCED VELOCITY= AXIAL INDUCED VELOCITY/VP 
< UT=NON- DIMENSIONAL TANGENTIAL COMPONENT OF 
INDUCED VELOCITY=TANGENT IAL INDUCED 
VELOCITY/VP 

GAMMA=NON-DIMENSIONAL Cl RCULATION=CIRCULATION/ 
(2»*PAI*VP*RP) 

=SQRT (1-Q**2) * SUMMATIONG (I) *U (1-1 r Q) , 

'I=T» 2, 3.,., 

REMARK: 1. ALL INPUT AND OUTPUT ARE IN NON-DIMENSIONAL 

FORMS 

2« 'INPUT DATA ARE MADE IN 

A. MAIN PROGRAM 

B. SUBROUTINE COMPLEX . 

C. SUBROUTINE JCNST1 

D. SUBROUTINE NOPTIML 

E. SUBROUTINE INDA12 

3. THIS PROGRAM CONSISTS OF A MAIN PROGRAM AND 
49 SUBPROGRAMS: 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1, • SUBROUTINE COMPLEX, 

2« SUBROUTINE MEANSQF, 3, SUBROUTINE JFDNC, 

4-, SUBROUTINE AECOEF, 5, ' SUBROUTINE JCNSTl, 

6r SUBROUTINE JCONSX., 7, SUBROUTINE JCEKl , 

8, SUBROUTINE JCENT, 9, SUBROUTINE AERO DIN, 

10c "SUBROUTINE NOPTIML, 11a, SUBROUTINE CTP12, . 
12. -SUBROUTINE AEROCF , 13. FUNCTION RAIJKM, 

14. SUBROUTINE CIRCU, 15, .SUBROUTINE UATCHBY, 
16. .SUBROUTINE HATIJ, 17. FUNCTION RAHDAF, 

18. . SUBROUTINE RAMCOE, 19, SUBROUTINE DOUCHB, 
20.. SUBROUTINE CHEBCF, 21. SUBROUTINE SUHODD, 
22. SUBROUTINE SUHCBB, 23. SUBROUTINE ZEROS, 
24, .SUBROUTINE DOUSUM,. 25. SUBROUTINE FACTOR,. 
26. -FUNCTION CIRCLF, 27. .SUBROUTINE INDVEL, 

28, SUBROUTINE VELOCIT, 2S. SUBROUTINE TCHEBY, 
30, SUBROUTINE OCHEBY, 31. - SUBROUTINE INDFACT, 
32, .SUBROUTINE INDA12, 33. FUNCTION GQDZ25, 

34. SUBROUTINE GQUDAD, 35* . SUBROUTINE IN DA 3, 
36* FUNCTION AIND3 , - 37* SUBROUTINE INDA4, 

38, FUNCTION AIND4, 39. SUBROUTINE COEAN, 

40. FUNCTION FAXIAL, 4! , : FUNCTION DLGAMA , 

42. SUBROUTINE WRENF, 43. SUBROUTINE MEANSQ, 
44, SUBROUTINE PRESUF, 45. SUBROUTINE FKGAT, 
46, SUBROUTINE RETARD, 47* SUBROUTINE PKO AT, 
48, FUNCTION FEHTON, 4 9, ’ SUBROUTINE ABSIAS. 


PROBLEM Is 

GIVEN: N BLADE, HUB 
VF 

RAMDAP ( OR FM AND TM) 

RAMDAI 

DETERMINE: AERODYNAMIC OPTIMUM CIRCULATION 


DISTRIBU TIOH 
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c 

c 

c 

c 

c 

c 

c 
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INDUCED VELOCITY COMPONENTS 
CT , CP, AND IDEAL EFFICIENCY 
ENSMSQ (OPTIONAL) 


PROBLEM II: 

GIVEN: NBLADE, HUB 
V? 

RAMDAP (OR FM AND TM) 

CT (OR CP) 

DETERMINE: AERODYNAMIC OPTIMUM CIRCULATION 
DISTRIBUTION 

INDUCED VELOCITY COMPONENTS 
CP (OR CT) AND IDEAL EFFICIENCY 
ENSMSQ (OPTIONAL) 

PROBLEM III: 

GIVEN: NBLADE, HUB 

VP 

RAMDAP (OR FM AND TH) 

CT (OR CP) 

THE TYPE OF CIRCULATION FUNCTION 

DETERMINE: NON-OPTIMUM CIRCULATION DISTRIBUTION 
INDUCED VELOCITY COMPONENTS 
CP (OR CT) AND IDEAL EFFICIENCY 
ENSMSQ (OPTIONAL) 


PROBLEM IV: 

GIVEN: NBLADE, HUB 

FM, TM, AND VF 
HUB 

UPPER BOUND AND LOWER BOUND OF CT 
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c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 
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r* 

c 

c 

c 

c 

c 

c 
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OPPEE BOUND AND LOWER BOUND OF CP 
UPPER BOUND AND LOWER BOUND OF RAH DAI 

DETERMINE: AEBOACOUSTIC OPTIMUM CIRCULATION 
DISTRIBUTION 

INDUCED VELOCITY COMPONENTS 
CT, CP, AND ETAI 
MINIMUM ENSEMBLE MEAN SQUARE OF 
ACOUSTIC PRESSURE 


USAGE: 

1. PROBLEM I: 

A. ENSMSQ IS NOT DESIRED 

INPUT: NBLADE, HOB 
VF, RAMDAP 
NTRY-1 

TRYRAM (1)=RAMDAI 

IPRBLM=1 

IMEAN=G 

B, ENSMSQ IS DESIRED 

INPUT: NBLADE, HUB 

VF, FH, TM (RAMDAP=FM/ (TM*VF) 
XDISTA, A ZIMUT {IN DEGREE) 
NTRY=1 

TRYRAM (1) =RAM DAI 

IPRBLM=1 

IMEAN=1 

2. PROBLEM II: 
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C 

c 

c 
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c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 
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c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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A. ENSMSQ IS SOT DESIRED' 

INPUT: NBLADE, HUB 

VP, RAMDAP (OR FH AND TH) 

CT (OR CP) 

IPRBLM=2 

IMEAN=Q 

NTRY=NUMBEE OF PRESET VALUES OF RAHDAI 

TRYRAM (I)=RAHDAI, I=1(1)NTRY 

CTP=CT IP CT IS SPECIFIED 

CTP=CP IF CP IS SPECIFIED 

NTP=1 IF CTP=CT 

N TP» NEj 1 IF CTP=CP 

B. ENSMSQ IS DESIRED 

INPUT: NBLADE, HUB 

VF , FH, TH 

XDISTA, AZIMUT (IN DEGREE) 

NTRY=NUMBER OF PRESET VALUES 
OF RAHDAI 

TRYRAH (I) = RAHDAI, 1=1 (1) NTRY 

CTP=CT IF CT IS SPECIFIED 

CT P=CP IF CP IS SPECIFIED 

NT P=1 IF CTP=CT 

NT P--> NE 0 1 IF CTP=CP 

IPRBLH=2 

IHEAN=1 

3, PROBLEH III: 

Ao -ENSMSQ IS NOT DESIRED 

INPUT: NBLADE, HUB 

VF, RAMDAP (OR FH AND TM) 

CTP=CT IF CT IS SPECIFIED 
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c 

r 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 
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CTP=CP IP CP IS SPECIFIED 

NTP=1 IF CTP=CT 

NTP « NE* 1 IF CTP = CP 

IPRBLM=3 

IHEAN=0 

GG{I) , 1=1 ,2,3,, . 

(SEE SUBROUTINE NOPTIML) 

WHERE 

GAHMA=K * SQR_T(1,-Q**2) * 

SUMMATION GG(I) *U (1-1 ,Q) ,, 1=1,2,,*. 
K: TO BE DETERMINED 

B. “ ENSMSQ IS DESIRED 

INPUT: NBLADS, HUB 
VF, FM, TM 
SDISTA, A2IMUT 
CTP=CP IF CT IS SPECIFIED 
CTP=CP IF CP IS SPECIFIED — 
NTP=1 IF CTP=CT 
NTP * NE* 1 IF CTP=CP 
IPRBLH=3 
IMSAN=1 

GG (I) , 1=1, 2, 3, .... 

(SEE SUBROUTINE NOPTIML) 

WHERE 

GAMHA=K * SQRT (1,-Q**2) * 

SUMMATION GG(I) *U (I-1,Q) 

1=1, 2,3, .;«••••* . 

K: TO BE DETERMINED 

4,: PROBLEM iv,: 



1 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

•c 

c 

c 

c 

c 


INPUT : NBL ADE, RUB 
VF, FM, TM 
XDIST A , AZIMUT 

UPPER BOUND AND LOWER BOUND OF CT 

UPPER BOUND AND LOWER BOUND OF CP 

UPPER BOUND AND LOWER BOUND OF RAMDAI 

AT. THE ZEROS OF T { JRAMDA+1 , Q) 

(FOR DETAIL, SEE SUBROUTINE COMPLEX) 
CTP=CT (OR CP) USED FOR FINDING 
THE FIRST APPROSIMATQN, 

(CT P MUST SATISFY CONSTRAIN) 

NTP=1 IF CTP=CT - - 

NTP . NE. 1 IF CTP=CP 

NTRY=NUMBER OF ASSUMED VALUES OF RAMDAI 

TRYRAM (I) =RAHDAI 1=1 (1) NTRY 

IPRBLM=4 

IMEAN=1 


PRECISION: SINGLE PRECISION 
REQUIRED DATA FILES 

COSGRA (SEE PROGRAM COECOS) NEEDED ONLY WHEN 

THE ALTERNATE METHOD DEVELOPED IN APPENDIX A 
IS APPLIED TO COMPUTE THE INDUCTION FACTORS 
PKOATD CONTAINS THE COEFFICIENTS, PKG , PKA, 

AND PKT, OF THE CHEBYSHEV EXPANSIONS OF 
THE KERNELS OF ACOUSTIC PRESURE 
NEEDED ONLY WHEN NTAPE . EQ. 2 
COMPUTED BY PROGRAM ITSELF IF NTAPE, NE. 2 

DESCRIPTION OF PARAMETERS 

HUB HUB RADIUS OF PROPELLER (I NPUT) 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NBLADE 

TM 

FH 

VF 


CT 

CP 

CTP 

ETAI 

JRAMDA 


MCI ECU 
ORAMDA (X) 
OQZERO(I) 


QCXRCU (X) 


XDISTA 

AZIMUT 

IDEG 


ICHEBY 


NUMBER CF PROPELLER BLADES (INPUT) 

TIP MUCH NUMBER OF PROPELLER (INPUT) 
FORWARD' MUCH NUMBER (INPUT) ' 
NON-DIMENSIONAL FORWARD VELOCITY OF 
PROPELLER Wo R* T REFERENCE VELOCITY 
VP, IF VP IS CHOSEN TO BE THE FORWARD 
VELOCITY OF PROPELLER, VF=1 (INPUT) 

THRUST COEFFICIENT 
POWER COEFFICIENT 
SPECIFIED CT OR CP (INPUT) 

IDEAL EFFICIENCY 

JEAMDA+1: NUMBER OF THE FUNCTION VALUES 
OF RAMDAI AT THE ZEROS OF T (JRAMDA+1 ,Q) 
THAT IS, THE NUMBER OF EXPLICIT 
INDENPENDENT VARIABLES (INPUT) 

NUMBER OF TERMS IN CHEBYSHEV EXPANSION OF 
THE CIRCULATION (INPUT) 

RAMDAI' (Q) AT ZEROS OF T (JRAMDA+1 ,Q) 

1=1 (1) JRAMDA+1 

ZEROS OF T (JRAMDA+1, Q) THAT IS, 


OQZERO (I)=COS ( (2*1+1) *PAI/ (2* (JRAMDA+1) ) 
ZEROS OF T (MCIRCU , Q) 

USED IN CALCULATION OF THE COEFFICIENTS 
OF THE CHEBYSHEV EXPANSION OF CIRCULATION 
NON-DIMENSIONAL DISTANCE BETWEEN THE 
PROPELLER AND THE OBSERVER (INPUT) 

AZIMUTH ANGLE OF THE OBSERVER IN DEGREE 
(INPUT) 

DEGREE OF THE DOUBLE CHEBYSHEV EXPANSION 
OF BOTH AXIAL AND TANGENTIAL INDUCTION 
FACTORS (INPUT) 

ICHEBY+1: NUMBER OF POINTS USED IN THE 

COMPUTATION OF THE COEFFICIENTS OF THE 
DOUBLE CHEBYSHEV EXPANSION FOR INDUCTION 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


QIND (I) 


EPSIND 

LUPIND 

KWENH 


NNOISE 

NOSCHB 


QNOISE (I) 
ETA (I) 

SQMEAN 

IMEAN 

IPHB1H 

NZERO 

NOS 

NTAPE 


ENS (I) 


FACTORS (INPUT) 

ZEROS OE T (ICHEBY+1 ,Q) , THAT IS 

QIND (I) =COS ( (2*1+1 ) *PAI/ (2* (ICHEBY+1) ) 
1=0 (1) ICHEBY 

LOWER LIMIT OF IA12 AND THE UPPER LIMIT 
OF I A3 (INPUT) 

2*PA I* LUPIND: UPPER LIMIT OF IA12 AND 
LOWER LIMIT OF IA4 (INPUT) 

CONTROL PARAMETER 

IF KWENH=1 , THE PROPELLER INDUCTION 
FACTORS ARE COMPUTED BY WRENCH'S 
FORMULA 

IF KWENH a NE. 1 THE PROPELLER INDUCTION 
FACTORS ARE COMPUTED BY ALTERNATE METHOD 
(INPUT) 

DEGREE OF THE CHEBYSHEV EXPANSIONS 
OF KO , KA, AND KT (INPUT) 

NOSCHB+1: NUMBER OF POINTS USED IN THE 
EVALUATION OF THE COEFFICIENTS OF THE 
CHEBYSHE V EXPANSIONS OF K0,KA, AND KT 
(INPUT) • 

ZEROS OF T (NOSCHB + 1, Q) 

ABSCISSAS OF THE 25-POINT GAUSS-LEGENDRE 
FORMULA 

MEAN SQUARE OF ACOUSTIC PRESSURE 
PROBLEM CONTROL PARAMETER (INPUT) 

PROBLEM CONTROL PARAMETER (INPUT) 

LENGTH OF ARRAY ZERO (INPUT) 

LENGTH OF ARRAY ORAKDA (INPUT) 

INPUT CONTROL PAREMETES (INPUT) 

NTAPE=0 IF DATA FILE PKOATD IS SOT 
PROVIDED 

NT APE=2 IF DATA FILE PKOATD IS PROVIDED 
BY USER 

ENS (1) =ENSMSQ CORRESPOND! NG TO 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


GT (I) 


GP (I) 


CTRAIN 

CPRAIN 


AERODYNAMIC OPTIMUM PROPELLER 
ENS (2) =ENSMSQ CORRESPONDING TO 
NON-OPTIMUM PROPELLER 


GT (1 )=COMPOTED CT 
AERODYNAMIC 
GT (2) =COKPUTED CT 
NON-OPTIMUM 
GP (1) =COHPOTED CP 
AERODYNAMIC 
GP (2) =COMPUTED CP 
NON-OPTIMUM 
MIN (GT (I) ) 

MAX (BP (I) ) 


CORRESPONDING 
OPTIMUM PROPELLER 
CORRESPONDING TO 
PROPELLER 
CORRESPONDING TO 
OPTIMUM PROPELLER 
CORRESPONDING TO 
PROPELLER 


DIMENSION ZERO (51) ,ORAHDA{15) , TRYRAM (5) 

COMMON /CON ST/PAI 

COMMON /DATAl/JRAMDA 

COMMON /DATA2/A0RAMD, ABAMDA (14) 

COMMON /DATA3/IDEG r ICHEB Y 
COMMON /DATA4/HA (3 1 , 31 ) , HT (3 1 , 31 ) 

COMMON /DATA5/QIND (51) 

COMMON /DATA6/HUB 

COMMON /DATA8/NBLADE,KRENH,LUPXHD 
COMMON /DAT A9/EP SIND 
COMMON /DATA1 Q/QCIRC U (15) 

COMMON /DAT A11/MCIRCU 
COMMON /DATA12/FM, TM, YF, RAMDAP 
COMMON /DATAl 4/G (1 5) 

COMMON /DATAl 5/CT # CP , ETAI 
COMMON /DATAl 6/OQZERO (15) 

COMMON /DATAl 7/XDI ST A/ AZIMUT, ANGLE, SANGLE, CANGLE 
COMMON /DATAl 9/QNOISE (51 ) 

COMMON /DATA20/NNOISE,NOSCHB 
COMMON /DATA21/ETA (25) 
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COMMON /DATA22/PK0 (25,31) ,?KA (25,31) ,PKT (25,31) 
COMMON /D AT A2 4/S QM EA N 
COMMON /DATA31 /I 73LDG,IVELCHB 
COMMON /DATA32/QVBL {51 ) 

COMMON /DATA33/UACHB (51) ,UTCHB(51) 

COMMON /DAT A34/CTRAIN, CP BA IN 
COMMON /DATA35/CX1 ,CT2 , CPI , CP2 
COMMON /DAT A36/I MEAN, IPRBLH 
COMMON /BATA37/EHS (2) 

COMMON /DATA38/XX (2, 15) 

COMMON /DATA39/GT (2) ,GP (2) 

PAID* 3 * 14159 2653589793 
C 

C*** INPUT data 

c 

c 

C PROPELLER DATA: 

C 

NBLADE=2 
HUB=0* 2 
C 

C AERODYNAMIC DATA: 

C 

vf*i. : 

FM=C. 2 
TM=0* 7426 
BAMBAF=FM/(TM* VF) 

C 

C ACOUSTIC DATA: 

C 

XDISTA=4*272 
AZIMUT=69* 44 
ANGLE=PAI*AZIMUT/1 80* . 

S ANGLE-SIN (ANGIE) 

CANGLE=COS (ANGLE) 


C 
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C PROBLEM CONTROL PARAMETER 

C 

IMEAN=1 

IPRBLM=4 

C 

C WRITE OUT THE PROPELLER DATA 

C 

60 FORMAT (///, 52 , H ******************************»') 

61 FORMAT (5X #'!***** *****»t) 

62 FORMAT (5X,"***** DEFINITIONS OF *****it) 

63 FORMAT (5X, "***** SYMBOLS «W*") 

64 FORMAT (51, »******************************»//) 

WHITE (6,60) 

WRITE (6,61) 

WRITE (6,62) 

WRITE (6 ,63) 

WRITE (6,61) 

WRITE (6,64) 

WRITE (6,65) v 

65 FORMAT (5X,"NBLADE= NUMBER OF PROPELLER BLADES"/,5X, 

$ "HUB=NQN-DIMENSIONAL HUB RADIUS'*/# 5X, 

$ "VF= NON- DIMENSIONAL ADVANCE SPEED OF PROPELLER *</#5X f 
$ »MF=ADVANCE MACH NUMBER OF PROPELLER"/# 5X, 

$ "MT=TIP MACH NUMBER OF PROPELLER"/, 5X,, 

$ "RAMDAP=MF/(MT*VF) ,REFERRENCS ADVANCE COEFFICIENT"/, 5X, 
$ "XDISTA=DISTANCE OF OBSERVER FROM THE PROPELLER"/ 

$ ,5X,"R=RADIAL COORDINATE OF PROPELLER BLADE SECTION"./# 

$ 5X,"Q= (2+R-HUB— -1) / (1-HUB) , -1 , SE. Q • LE« 1»/,5X, 

$ "CT=THRUST COEFFICIENT"/# 5X, 

$ »CP=POWER COEFFICIENT”/, 5X, 

$ "G (I) COEFFICIENTS OF THE CHEBYSHEV EXPANSION"/, 52, 

$ "OF CIRCULATION, I. E. GAMMA"/, 5X, 

$ "UA=AXIAL COMPONENT OF INDUCED VELOCITY"/, 5X, 

$ «UT=TANGENTIAL COMPONENT OF INDUCED VELOCITY"/, 5X, 

$ "RAMDAI=HYDRODYNA MIC ADVANCE COEFFICIENT"/, 5X, 

$ »GAMHA=NON-DIMENSIONAL CIRCULATION"///) 
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50 FORMAT {/// ,5X t **************************************** ”) 

51 FORMAT (5X, "***** *****'') 

52 FORMAT (5X ,«***** PROPELLER DATA AND *****»') - 

53 FORMAT (5X, « ***** PROGRAM CONTROL PARAMETERS *****"} . 

54 FORMAT (5X, m*************************************** *'//) 
WRITE (6 #50} 

WRITE (6,51) • 

WRITE (6,52) 

WRITE (6,53) 

WRITE (6,51) 

WRITE (6,54) 

WRITE (6, 19} NBLADE,HBB,VF,FM,TH,RAMDAP, 

$ XDISTA ,AZIMUT 

19 FO RM AT (//, 5 X , » N BLADE-*’ ,12/ , 5X , »HUB=» ,F1 0* 5/5 X, 

$ »VF= » , F 1 0 . 5/5 X , " MF=» , FI 0. 5/5X , »H T= H , F 1 0. 5/5X , 

$ »EAMDAP=" , F10 • 5/, 5X , « XDISTA=" ,F1 0 * 5/5X, 

$ “AZIMUTH ANGLE=", FIG, 5, 2X,” (DEGREE) "//)- 
C ’ 

C 

C PARAMETERS USED IN THE COMPUTATION OF 

C INDUCTION FACTOR 

C 

KWENH=1 
EPSIN D=0o 
LUPIND=4 
C 
C 

C PARAMETERS FOR NUMESCIAL SOLUTION OF 

C OPTIMAL CIRCULATION 

C 
C 

TRYRAM (1 ) =0* 315 
TRYRAM (2) =0.32 
TRYRAM (3) =0.325 
NTRY=3 
NTP=1 
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CTP=0» 3556 

NRAM=5 

URAMDA=9 

IDEG=10 

ICHEB Y=10 

HCIRCU=10 

NNOISE=10 

NOSCHB=10 

IVELDG=2*IDSG+MCIRCU-1 

IVELCHB=IVELDG 

C 

C PRINT OUT THE PROGRAM CONTROL PARAMETERS 

C 

WRITE (6 ,20} JRAHDA, IDEG,ICHEBY , MCIRCU , 

$ NNOISE , NOSCHB 

20 FORMAT (5X,»JRAMDA=»,I3/5X,"IDEG=",I3/5X, 

$ »ICHEBY= ,, ,I3/5X, ,, MCIRC0= ,, / I3/5X, 

$ "NNOISE =»,I3/5X,"N0SCHB=»,I3//) 

C 

c 

C PROGRAM PARAMETERS 

C 

NTAPE=2 

NZERO—51 

NOR=15 

C 

C*** END OF INPUT 
C 

C , COMPUTE QIND-(I) THE ZEROS OF T (ICHEBY+1 , Q) 
C 

NDEGHE= ICHEBY+1 

CALL ZEROS (NDEGRE, NZERO , ZERO) 

DO 1 1=1, NDEGRE 
QIND (I) = ZERO (I) 

1 CONTINUE 
C 
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C COMPUTE QCIECU(I), THE ZEROS OF T(MCIRCU,Q) 

C 

NDEGRE=MCIRCU 

CALL . ZEROS (NDEGRE, NZERO, ZERO) 

DO 2 1=1, NDEGRE 
QCIRCU (I)=ZERO (I) 

2 CONTINUE 
C 

C COMPUTE OQZERO(I), THE ZEROS OF T (JRAMDA+1 ,Q) 

NDEGRE=JRAMDA+1 
CALL ZEROS (NDEGRE, NZERO, ZERO) 

DO 3 1=1, NDEGRE 
OQZEEO (I) =ZERO (I) 

3 CONTINUE 
C 

C COMPUTE THE ZEROS OF T (NOSCHB+1 ,Q) 

C 

NDEGRE=NOSC HB+1 

CALL ZEROS (N DEGRE, NZERO , ZERO) 

DO 4 1=1 , NDEGRE 
QNOISE(I) =ZERO (I) 

4 CONTINUE 
C 

C COMPUTE THE ZEROS. OF T (I VELCHB+1 ,Q) 

C 

NDEGRE=IVELCHB+ 1 

CALL ZEROS (NDEGRE, NZERO, ZERO) 

DO 12 1=1, NDEGRE 
QVEL (I) =ZERO (I) 

12 CONTINUE 
C 

IF (IMEAN . NE» 1} GO- TO 5 
C 

C CALL SUBROUTINE ABSIAS TO OBTAIN ETA (I) 


CALL ABSIAS 
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C 

C CALX. SUBROUTINE PK0AT 10 OBTAIN PKQ (I, J) , 

C PKA {I , J) , AND PKT (X/ J) t 1=1(1)25, J=1 (1) NNOISE+1 

C IF NTAPE=0 

C 

IF (NTAPE.NE.O) GO TO 3 
• CALL PKGAT 
GO TO 5 
C 

C READ IN PKG (I, J) , PKA (I,J) , AND PKT (I,J) 

C PROS DATA FILE PK0ATD 

C 

8 REWIND2 
JEND=NN0ISE+1 
DO 9 1=1 ,25 

DO 10 J=1,JSND 

READ (2,11) PKO (I,J) ,PKA (I,J) ,PKT (I,J) 

.10 CONTINUE 

9 CONTINUE 

11 FORMAT (3 (2X,E23« 16)) 

C 

C CALL SUBROUTINE AERO DIN TO OBTAIN THE 

C INITIAL SET OF POINTS USED IN THE 

C SUBROUTINE COMPLEX 

C 

c 

c 

5 IF {IPR8LM. HE, 3) 

$ CALL ABRQDYN (TRYRAM ,NRAM,ORAH DA , NOR , NTKX, NT P,CT P) 
IPR« {IPRBLM-1 ) * (IPRBLM-2) 

IF(IPB.NE*0) CALL NOPTIML (ORAMDA , NOR, NT P ,CTP) 

IF (IPRBLM , NB. 4) STOP 
C 

C CALL COMPLEX TO OBTAIN THE AEROACOUSTIC 

C OPTIMAL CIRCULATION DISTRIBUTION 

C 
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CALL COMPL EX'(ORAMDA, NO R) 

CALL RAMCOE (OR AMDA ,NOR) 

CALL HATIJ 
CALL CIRCUL 
CALL UATC BBY 

40 FORMAT (//, 5X, »*********************************") 

41 FORMAT (5X, ****** *****»«) 

42 FORMAT (5X, »**♦** ASROACOUSTIC SOLUTION *****>') 

43 FORMAT (5X, «*********************************»//) 

WRITE (6,40) 

WRITE (6, 41). 

WRITE (6,42) 

WRITE (6,41) 

WRITE (6,43) 

WRITE (6 ,31) HCIRCU 

31 FORMAT (5X, ’’NON- DIMENSIONAL CIRCULATION”//, 5X , 

$ "GAMMA (Q) = SQRT (1“Q**2) +SUHMATION G (I) * 0 (1-1 ,Q) , 1=1," 
$ , 12//, 5X,» WHERE G(I):»//) 

WRITE (6,32) (I,G (I) ,I=1,MCIRCU) 

3 2 FORMAT (5X, "G (»,I2, M ) =",F25,10) 

CALL AEROCF 

WRITE (6,33) CT,CP,ETAI 

33 FORMAT (//,5X, "CT=» ,F10 * 5//, 5X, »CP=» , F10, 5//,5X, 

$ "IDEAL EFFICIENCY*" ,F 10. 5//) 

WRITE (6,44) • 

44 FORMAT (///, 1 OX, » ++R++" , 1 2X , " ++Q++ " ,1 1 X, "++UA++" , 

$ 11X, "++UT++ " ,9X ,"++ RAMDAI++" , 8X ,' "++GAMM A ++"//} 

R=HUB 

3 5 IF (R* GT. 1. ) GO TO 37 

Q~ (2. *R- (1* +HUB) ) / (1 • - HUB) 

CALL 7EL0CIT (UA,UT,Q) 

RAH=R* (VF + UA)/ (R/RAMDAP+UT) 

WRITE (6,36) R,Q,UA,UT,RAM,CIHLF(Q) 

R=R+Q, 1 
GO TO 35 

FORMAT (6(2X,F15«7) ) 


36 



STOP 

END 


*****************2**:*** *************** 
***** ***** 

***** COMPLEX METHOD ***** 
***** ***** 
******************** ***************** 


SUBROUTINE COMPLEX (ORAHDA, NOE) 

PURPOSE 

TO FIND THE AEROACOUSTIC OPTIMAL 
CIRCULATION 

DESCRIPTION OF PARAMETERS 

SEE SUBROUTINE JCONSX AND MAIN PROGRAM 


DIMENSION ORAMDA (NOR) 

DIMENSION X (1 5 ,1 5) ,R (15,15) ,F (15) 
DIMENSION G (15) ,H (15) ,XC (15) 
COMMON /DATA12/FH, TM,VF, RAMDAP 
COMMON /DAT Al 6/OZERO (1 5) 

COMMON /JTRAIN/JYES 
COMMON /DATAl/JRAMDA 
COMMON /DATA6/HUB 
COMMON /DATA15/CT, CP, ETAI 
COMMON /DATA30/ENSMSQ (15) 

COMMON /DATA34/CTRAIN, CPRAIN 
COMMON /DATA37/ENS (2) 
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COBMON /DATA38/XX(2, 15) 

COMMON /DATA39/GT (2) ,GP (2) 

INTEGER GAMMA 

C ' 

C*** input PARAMETERS 
C 

KO=6 
KX=15 
LX=1 5 
HX=15 
NX=1 5 

N=J RAMDA+1 
M=N+1 
L=fl 
K=N+2 
ALPHA=1 a 3 
BETA=1*E-8 
GAMMA=4 
DEITA=0.0001 
ITH AX=1 40 
C 

C *** END OF INPUT 
C 

c 

c 

21 FORMAT (//,5 X fH**'*************'*********'**********) 

22 FORMAT (5X , ******* . *****«») 

23 FORMAT (SX,. 11 ***** BEGIN COMPLEX METHOD *****««) 

24 FORMAT (5X, «'********************************»//) 
WRITE (6,21). 

WRITE (6,22) 

WRITE (6,23) 

WRITE (6 ,22) 

WRITE (6,24) 

C 

C INPUT THE INITIAL SET OF POINTS 
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DO 1 1=1 ,N', 

X (1,I)=XX(1,I) 

X (2,1) =XX (2,1) 

1 CONTINUE 

X (1 ,N + 1 ) =GT (1) 

X (1,N+2)=GP (1) 

X (2,N + 1)=GT (2) 

X (2,N+2)=GP (2) 

C 

C CONSTRAINT ON CT OH CP 

C 

IF (GT (1)*LE,GI (2)J CTEAIN=GT(1) 

IF (GT (1) , GE, GT (2) ) CIHAIN=GT (2) 

IF (GP (1) .IE. GP (2) ) CPHAIN=GP(2) 

IF (GP(I).GE, GP (2) ) CPRAIN=GP(1) 

C 

C CHECK THE INITIAL SET OF POINTS 

C 

DO 40 IW=1,2 

IU=IW 

JYES=1 

CALL JCNST1 (N,M,K,X, G,H,IU,L,KX,LX,aX,NX) 

DO 2 1=1, B 

IF ( (X (IW ,1) -G (I) ) o LT* 0* ) GO TO 3 
IF ( (H (I)-X(IW,I) )„ LT.O.) GO TO 3 

2 CONTINUE 
40 CONTINUE 

GO TO 4 

3 WHITS (6,5) 

5 FORMAT (//,5X,»*** EXECUTION IS TERMINATED**"//) 
WHITE (6,6) 

6 FORMAT (5X,"THE INITIAL SET OF POINTS THAT»A5X, 
$ "DOESE NOT SATISFIES CONSTRAINTS"//) 

DO 8 1=1, M 

WRITE (6 , 9) I W,I ,X (1 ,1) 



8 CONTINUE 

9 FORMAT (5X,»X(",I2, 1, , ,, ,I2, M )="»F15.7) 

STOP 

C 

C GIVE ENSMSQ(I) AND ENSMSQ(2) 

C 

4 ENSMSQ (1) =ENS (1) 

ENSMSQ (2) =ENS (2) 

C 

C INPUT THE RANDOH NUMBERS R (K,N) 

C 

C UES SUBPROGRAM FUNCTION RANF (A) TO GENERATE 

C RANDOH NUMBERS. 

C 

DO 30 1=1, K 
DO 31 J=1,N 
R (I ,J)=RANF (AA) 

31 CONTINUE 

30 CONTINUE 

C 

C CALL SUBROUTINE JCONSX TO OBTAIN THE 

C AEROACOUSTIC SOLUTION 

C 

CALL JCONSX (N, H, K, IT MAX, ALPHA, BETA , GAMMA , 

$ DELTA ,X , R , F, IT , IEV2 ,KO,G,H,XC,L, KX,LX , MX, NX) 

DO 10 1=1, N 
ORAMDA (I) =X (IEV2,I) 

10 CONTINUE 

CT=X (IEV2,N+1) 

CP=X(IEV2,N+2) 

. SQMEAN=-F (IEV2) 

ETAI=VF*CT/CP 
WRITE (6,11) 

1 1 FORMAT (//, 51, H**#*******#**************#******^ 

12 FORMAT (5X," ***** *****iij 

13 FORMAT (5X, ”***** OUTPUT FROM COMPLEX *****'•) 
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1 8 FORMAT {5X , " ***********#*************+*******) 

WRITE (6,12) 

WRITE (6,13) 

WRITE (6,12) 

WRITE (6,18) 

IDG=JRAMDA+1 
WRITE (6,14) 

14 FORMAT {//,5X, "HYDRODYNAMIC ADVANCE COEFFICIENTS"//) 

WRITE (6,20) 

20 FORMAT (//, 1 5X, "++R++" , 12 X,»++Q++» ,9 X, "++RAMDAI++"//) 

DO 15 I— 1’,IDG- 

ROU=0 1> 5* ( (1 .-HOB) *OZERO (I) +1. +HOB) 

WRITE (6,16) ROD ,OZERQ (I) ,ORAMDA ( I) 

15 CONTINUE 

16 FORMAT (5X,3 (2X,F15,7) ) 

NRITE-{6 , 1 7) CT , C P , ET A-I ,-S QUEAN 

17 FORMAT (//,52, ,, CT=" ,F15. 7//,5X, "CP=",F15*7//,5X, 

$ "IDEAL EFFICIENCY^ 1 , F15*7//, 5X, "ENSEMBLE MEAN SQUARE=», 
$ E25.12//) 

RETURN 

END 

C 

n 

*«*» 

SUBROUTINE MEANSQF (X , I ,11 , KX , LX) 

C 

C PURPOSE 

C 

C TO COMPUTE THE OBJECTIVE FUNCTION 

C BY CALLING SUBROUTINE MEMSQ AND STORE 

C THE RESULT IN ARRAY ENSMSQ(I) FOR THE 

C ITH SET OF POINTS 

C 
C 

DIMENSION X(KX,LX) 

COMMON /DATA3G/ENSMSQ ( 15) 

COMMON /DATA2 4/SQMEAN 
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CALL MEANSQ 
ENSMSQ {I)=SQHEAN 
RETURN 
EKD 
C 
C 

SUBROUTINE JFUNC (N ,M , K , X , F r I , L ,KX, LX, MX, NX) 

C 

C 

C PURPOSE 

C 

C 

C TO SUPPLY THE FUNCTION AL VALUE 

C OF THE OBJECTIVE FUNCTION THAT IS 

C TO BE MAXIMIZED 

c - - - 

c 

DIMENSION X (KX/LX) ,F (KX) 

COMMON /DATA3G/ENS MS Q { 15) 

P (I) =-ENSMSQ(I) 

RETURN' 

END 

C 

C‘ 

SUBBOUTINE AECOEF (N r M, K,X,I,KX,LX) 

C 

C PURPOSE 

C 

C TO COMPUTE CT AND CP FOR A GIVEN SET OF POINTS 

C 

DIMENSION ORAMDA (15) ,X (KX,LX) 

COMMON/DATA1 /JRAMDA 
N0R=15 
JX=JRAHDA+1 
DO 1 J=1 r JX 
ORAMDA (J)=X (I,J) 
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1 CONTINUE 

CALI RAHCOE (ORAMDA,NOE) 
call' HATIJ 
CALL CIRCUL 
CALL U AT CHS I 
CALL ABROCF 
BETUSN 
END 
C 
C 

SUBROUTINE JCNST1 (N, H , K ,X , G, H,I, L, KX, 

$ LX, ax, NX) 
c 
c 

C PURPOSE 

C 

C 

C TO PROVIDE THE CONSTRAINTS 

C 

C 

DIMENSION X (KX,LX) ,G (MX) ,H (MX) 

COHHON /BATA15/CT,CP,STAI 
■COHHON /JTRAIN/JYES 
. COHHON /DAT A3 4/CTRAIN, CPEAIN 
C 
C 

C*** INPUT THE CONSTRAINTS 
C 

C EXPLICIT CONSTRAINTS (CONSTRAINTS ON RAHDAI) 

C 

DO 1 ITAB— 1 , N 
G (ITAB) =0, 1. 

H (ITAB) =0,6 
1 CONTINUE 
C 

C*** END OP INPUT' OF EXPLICIT CONSTRAINTS 
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c 

IF (JIES . HE. 1) RETURN 
C 

C ■ TO PROVIDE THE VALDES AND CONSTRAINTS 
C OF THE IMPLICIT VARIABLES 

C 

c 

C CALL AECOSP TO OBTAIN CT AND CF 

C 

CALL AECOEF (N,H,K,X,I, KX,LX) 

X (I,N+1)=CI 
X {I, N+2) -CP 
C 

C*** INPUT THE IMPLICIT CONSTRAINTS 
C (CONSTRAINTS ON CT AND CP) 

C 

S (N+1) =CTSAIN-1. 1-6 
H (N+1) =CT+1 * ' 

G (N+2) =0,5 
H (N+2) =5. 

C 

C*** END OF INPUT OF IMPLICIT CONSTRAINTS 

C 

C 

c***' END OF INPUT OF CONSTRAINTS 
C 

RETURN 

END 

C 

C 

C 

C 

SUBROUTINE JCQNSX (N, II, K, IT MAX, ALPHA, BETA , GAMMA, 
S ‘DELTA, X,R,F, IT, IEV2,K0,G, H,XC,L, KX,LX, MX,NX) 

C 

c 
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C THIS IS A MODIFIED PROGRAM OF RICHARDSON AND 

C KUESTEES* ORIGINAL PROGRAM (COMM, ACM 16 , 

C ' AUG, L97 3, PP, 487-489) 

C THIS MODIFIED PROGRAM OSES TWO INITIAL 

C FEASIBLE SOLOTIONS (AERODYNAMIC OPTIMUM SOLUTION 

C AND NON-OPTIMUM SOLUTION) AS INPUT AND' 

C COMPUTE THE IMPLICIT VARIABLES ONLY WHEN THEY ARE 

C NEEDED 

C 

C PURPOSE 

C 

C TO FIND THE CONSTRAINED MAXIMUM OF A FUNCTION OF 

C SEVERAL VARIABLES BY THE COMPLEX METHOD OF H. 'J.-BOX*: 

C THIS PROGRAM IS WRITTEN BY JOEL A, RICHARDSON AND 

C J, -L, KU ESTER, COMM, . ACM 16 (AUG, 1 973) , (487-489) 

C THIS IS THE PRIMARY SUBROUTINE AND COORDINATES THE 

C SPECIAL PURPOSE SUBROUTINES (JCEKl , JCENT, 

C JFUNC, JCNSTl) o .INITIAL GUESSES OF THE INDEPENDENT 

C VARIABLES, RANDAOM NUMBERS, 

C SOLUTION' PARAMETERS, DIMENSION LIMITS AND PRINTER CODE 

C DESTINATION ARE OBTAINED FROM THE MAIN PROGRAM 

C (SUBROUTINE COMPLEX) . 

C FINAL FUNCTION AND INDEPENDENT VARIABLE VALUES ARE 

C TRANSFERRED TO THE MAIN PROGRAM FOR PRINTOUT, . 

C INTERHINATE PRINTOURS ARE PROVIDED IN THIS SUBROUTINE. 

C THE USER MUST PROVIDE THE MAIN PROGRAM AND THE 

C SUBROUTINES THAT SPECIFY THE FUNCTION (JFUNC) 

C AND CONSTRAINTS (JCNSTl) , 

C FORMAT CHANGES 'HAY BE REQUIRED WITHIN THIS SUBROUTINE 

C DEPENDING ON THE PARTICULAR PROBLEM UNDER CONSIDERATION. 

C USAGE 

C CALL JCONSX (N,M,K,ITAX, ALPHA, BETA, GAMMA,DELTA,X,R,F, IT 

C ,IEV2,KO,G,B,XC,L,KX,LX,HX,NX) 

C SUBROUTINE REQUIRED 

C JCEKl (N,M,K,X,G,H,I,KODE, XC, DELTA , L , Kl , KX,LX ,NX ,NX) 

C CHECKS ALL POINTS AGAINST EXPLICIT AND IMPLICIT 
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C 

C 

c 

c 

c 

c 

c 

c 

.c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

€ 


CONSTRAINTS AND APPLIES CORRECTION IF VIOLATIONS ARE 
FOUND 

JCEHT (N,M,K,IEV1.,I,XC, X,L,K1 ,RX,.LX,MX,NX) 

CALCULATES THE CENTROID OF POINTS 
JFUNC (N,H,K,X,F,I,L,KX,LX,MX,NX) 

SPECIFIES THE OBJECTIVE FUNCTION (USER' SUPPLIED) 

JCNSTl (N,M,K,X,G,H,I # L,KX,LX,MX,NX) 

SPECIFIES EXPLICIT AND IMPLICIT CONSTRAINT LIMITS 
(USER SUPPLIED), ORDER EXPLICIT CONSTRAINTS FIRST 
DESCRIPTION- OF THE PARAMETERS 

KX NUMBER OF ROSS. IN THE DIMENSION STATEMENT 

FOR X, F, AND R , KX« GE* K-DEFINE IN MAIN PROGRAM 
(USER SUPPLIED) 

LX NUMBER OF COLUMNS IN THE DIMENSION STATEMENT 
FOR X, LX. GE» L- DEFINE IN MAIN PROGRAM (USER 
SUPPLIED) 

MX DIMENSION OP G AND H, MX, GE. M-DEFINE IN MAIN 
PRQGBHA (USER SUPPLIED) 

NX DIMENSION OF‘ XC AND THE NUMBER OF COLUMNS 

NX, GE.N-DSFINE IN MAIN PROGRAM {USER SUPPLIED} 

N NUMBER OF EXPLICIT INDEPENDENT VARIABLES- 
DEFINED' IN THE MAIN PROGRAM 
M NUMBER OF SETS OF CONSTRAINTS -DEFINED IN 
THE MAIN PROGRAM 

K NUMBER OF POINTS IN THE COMPLEX -DEFINED IN 
THE MAIN PROGRAM 

ITMAX MAXIMUM NUMBER OF ITERATIONS -DEFINED IN MAIN 
PROGRAM 

ALPHA REFLECTION FACTOR - DEFINED IN MAIN PROGRAM 
BETA CONVERGENCE PARAMETER -DEFINED IN MAIN PROGRAM 
GAMMA CONVERGENCE PARAMETER - DEFINED IN MAIN PROGRAM- 
DELTA EXPLICIT CONSTRAINT VIOLATION 
-DEFINED IN MAIN PROGRAM 

X INDEPENDENT VARIABLES - DEFINED IN MAIN PROGRAM 
H RANDOM NUMBERS BETWEEN 0 AND 1 - DEFINED IN MAIN 
PROGRAM 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

n 

c 

c 

c 

c 

c 

c 


10 

20 


F OBJECTIVE FUNCTION ** DEFINED IN SUBROUTINE JFUNC 

IT ITERATION INDEX -DEFINED IN SUBROUTINE JCONSX 

ISV2 INDEX OF POINT WITH MAXIMUM FUNCTION VAIUV 
-DEFINED IN SUBROUTINE JCONSX 
IEV1 INDEX OF POINT WITH HINTS UH FUNCTION VALUE 
-DEFINED IN SUBROUTINE JCONSX 
KO PRINTER UNIT NUMBER -DEFINED IN MAIN PROGRAM 
G LOWER CONSTRAINT -DEFINED IN SUBROUTINE JCNST1 

H UPPER CONSTRAINT -DEFINED IN SUBROUTINE JCNST1 

XC CENTROID -DEFINED IN SUBROUTINE JCENT 
L TOTAL NUMBER OF INDEPENDENT VARIABLES {EXPLICIT + 
IMPLICIT) -DEFINED IN MAIN PROGRAM 
I POINT INDEX -DEFINED IN SUBROUTINE JCONSX 
NODE KEY USED TO DETERMINE IF IMPLICIT CONSTRAINTS 
ARE PROVIDED -DEFINED IN SUBROUTINE 
JCONSX AND JCEKl 

K1 DO LOOP LIMIT -DEFINED IN SUBROUTINE JCONSX 
JYES IF. JYES=1, IMPLICIT VARIABLES ARE COMPUTED 
IN- SUBROUTINE JCNST1 

IF JYES • NE. '1 IMPLICIT VARIABLES ARE NOT COMPUTED 


DIMENSION X (SX,LX) , R {KX, NX) ,F (KX) ,G{MX) ,H (MX) ,XC (NX) 

INTEGER GAMMA 

COMMON /JTRAIN/JYES 

JYBS=2 

IT=1 

WRITE {KO, 99995) IT 
KODE=0 

IF {H-N) 20,20, 10 

K0DE=1 

CONTINUE 

DO 40 11=3, K 

DO 30 J=1,N 

X{II,J) =0. 
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30 CONTINUE 
40 CONTINUE 

C CALCULATE COMPLEX POINTS AND CHECK AGAINST CONSTRAINTS 

DO 60 11=3 , K 
DO 50 J=1 , N 
I=II 

CALL JCNSTl {N,M,K,X, G,H,I,L, KX,LX,HX,NX) 

X (II,J)=G(J)+R (II,J)*(H(J) **G (J) ) 

5C CONTINUE 
K1 = II 

CALL JCEKl (N,M,K,X,G,H,I, KOBE, XC, DELTA ,L, Kl , 

$ KX, LX, MX, NX) 

WRITE (KO, 99999) II, (X (II, J) ,J=1 ,N) 

60 CONTINUE 
Kl=K 

DO 70 1=1, K 

CALL JFUNC (N,M,K,X,F,I,L,KX,LX,MX,NX) 

70 CONTINUE 
KOUNT=1 
IA=0 

C ' FIND POINT WITH LOWEST FUNCTION VALUE 
WRITE (KO, 99998) (F(I),I=1,K) 

WRITE (KO, 99993) (X(I,N+1) ,1 = 1, K) 

WRITE (KO, 99994) (X(I,N+2) ,1=1, K) 

80 IEV1 = 1 

DO 100 ICM=2,K 

IE (F (IEV1) -F (ICM) ) 100,100,90 
90 IEV1=ICM 
100 CONTINUE 

C FIND POINT WITH HIGHTEST FUNCTION VALUE 

IEV2=1 

DO 120 ICM=2,K 

IF (F(IEV2) -F (ICS)) 110,110,120 
110 IEV2=ICM 
120 CONTINUE 

C CHECK CONVERGENCE CRITERIA 
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IF (F (IEV2) - (F (IEVl ) +BETA) ) 140,130,130 
130 KOUNT=1 

GO TO 150 
140 KOUNT=KOUNT+1 

IF (KOUNT-GAMMA) 150,240,240 
C REPLACE POINT BY LOWEST FUNCTION VALUE 

15C CALL JCENT (N,M,K,IEVl,I, XC,X,L,K1 ,KX,LX,MX,NX) 

DO 160 J=1,N 

X (IEV1 ,J) = (1.+ALPHA) *(XC (J) ) -ALPHA* (X (IEVl,J)) 

160 CONTINUE 
I=IEV 1 

CALL JCEKl (N, H,K,X,G,H,I, KOBE, XC, DELTA, L,K1, 

$ KX, LX, MX, NX) 

CALL JFUNC (N,H,K,X,F,I,L,KX,LX,HX,NX) 

C REPLACE NEW POINT IF IT REPEATS AS LOWEST FUNCTION VALUE 

170 IEV2=1 

DO 190 ICN=2,K 

IF (F (IEV2) -F (IC£t) ) 190,190,180 
180 IEV2= ICM 
190 CONTINUE 

IF (IEV2-IEV1) 220 ,200 ,220 
200 DO 210 ' JJ-1 ,N 

X (IEVl , JJ) = (X (IEV1 ,JJ) +XC (JJ) ) /2, 

210 CONTINUE 
I=IEV1 

CALL JCEKl {N,M,K,X,G,H,I, KODE, XC, DELTA ,L, Kl , 

$ KX, LX, NX, NX) 

CALL JFUNC (N,M,K,X,F,I,L,KX,LX,KX,NX) 

GO TO 170 
220 CONTINUE 

WRITE (KO, 99997) {X (IEVl, JB) ,JB=1,N) 

WRITE (KO, 99998) (F(I),I = 1,K) 

WRITE (KO, 99993) (X (1,11+1) ,1=1 ,K) 

WRITE (KO, 99994) (X(I,N+2) ,1=1, K) 

WRITE (KO, 99996) (XC (J) , J=1,R) 

IT=IT+1 
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IF(IT-ITMAX) 230,230,240 
230 CONTINUE 

WRITE (KO, 99995) IT 
GO TO 80 
240 RETURN 

99999 FORMAT (1H , 15X, 21H COORDINATES AT POINT, 

$ 14/5 (FI 5* 7 , 2X) ) 

99998 FORMAT (1H , 2GX, 16H FUNCTION VALUES, /5 (El 5, 7 , 2X) ) 
99997 FORMAT (IE , 20X,16H CORRECTED POINT, /5 (F15.7,2X) ) 

99996 FORMAT (IE , 21H CENTROID COORDINATES, 2 X, 5 (FI 5. 7,2X) } 
99995 F0EMAT.(1H , //I OH ITERATION ,4 X, 15) 

99994 FORMAT (1 H , .2DX,- "POWER COEFFICIENTS", /5 (F15. 7,2X) ) 
99993 FORMAT (IE , 2DX,- "THRUST COEFFICIENTS", /5 (F15 a 7, 2X) ) 

END 

C 

c 

SUBROUTINE JCEKl (N ,M ,K ,X ,G,H, I,KODE,XC , DELTA , L, K1 , 

C 

C 

$ KX, IX, MX, NX) 

C PURPOSE 

C TO CHECK ALL POINTS AGAINST THE EXPLICIT AND IMPLICIT 

C CONSTRAINTS AND TO APPLY CORRECTIONS IF VIOLATIONS 

C ARE FOUND 

C 
C 

C USAGE 

C 

C 

C CALL JCEKl (N,M,K,X,G, H, I, KOBE, XC, DELTA, L, Kl , KX,LX, MX, NX) 

C SUBROUTINE REQUIRED 

C JCENT (N,M,K,IEV1 ,I,XC,L,Kl ,KX,LX,MX,NX) 

C JCNST1 (N,M,K,X,G,H,I,L,KX,LX, MX, NX) 

C DESCRIPTION OF PARAMETERS 

C PREVIOUSLY DEFINED IN SUBROUTINE JCONSX 

DIMENSION X(KX,LX) ,G (MX) ,H (MX) ,XC (NX) 
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COMMON /JT RAIN/JYE S 
10 KT=0 
JYES=2 

CALL JCNST1 (N,M,K,X,G,H,I, L, KZ,LX,MX,NX) 

C CHECK AGAINST EXPLICIT CONSTRAINTS 

DO 50 J=1,N 

IF (X (I,J)-G (J) ) 20 ,20, 30 
20 X (I, J) =G (J) +DE1TA 
GO TO 50 

30 IF {H (J) -X {I , J) ) 40,40,50 
40 X (I,J)=H (J) -DELTA 
50 CONTINUE 

IF (KODE) 110,110,60 

C CHECK AGAIN STTHE IMPLICIT CONSTRAINTS 

60 CONTINUE. 

NN=N+1 

JXES=1 

DO 100 J=NN,M 

CALL JCNSTl (N,H,K,X,G,H,I,L,K2 ,LX , MX , NX) 

IF (X (I , J) -G (J) ) 80,70,70 
70 IF (H (J)-X (I,J) ) 80,120,120 
80 IEVl=I 
KT=1 

CALL JCENT (N,M,K,IEV 1 ,I,XC ,X,L,Kl ,KX,LX, MX, NX) 
DO 90 JJ=1,N 

X(I,JJ)= (X (I,JU)+XC(JJ) )/2* 

90 CONTINUE 
JYES—1 
GO TO 100 
120 JYES=2 
100 CONTINUE 

IF (KT) 110,110,10 
110 CALL MEANSQF (X ,1 ,N ,KX ,LX) 

JYES=2 

RETURN 

END 



SUBROUTINE JCENT (N ,H ,K , IEV1 ,1 , XC, X,L, Kl ,KX ,LX,MX ,NX) 


PURPOSE 

TO CALCULATE THE CENTROID OF POINTS 
USAGE 

CALL JCENT (N,H , K , IEVl , I , XC , 2, L,Kl , KX,LX,HX,NX) 

SUBROUTINE REQUIRED 

NONE 

DIMENSION X (KX /LX) ,XC (NX) 

DO 20 J=1 , N 
XC(J)=0, 

DO TO IL=1 , Kl 
XC (J)=XC (J) +X [IL , J) 

CONTINUE 

RK=K1 

XC (J) = (XC (J)-X (IEV1,J) )/ (RK-1*) 

CONTINUE 

RETURN 

END 


** jjc * «^k if * if jjt jp ^Ck^k «^k^C «jc 

***** AERODYNAMIC MODEL ***** 
***** ***** 
******************** **** ***** 


SUBROUTINE AERODYN (TRY RAM, NRAH,ORAMDA, 
$ NOR, NTRY,NTP,CTP) 


PURPOSE 
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C 

C TO ESTIMATE THE FIRST APPROXIMATION OF THE 

C HYDRODYNAMIC ADVANE COEFFICIENT BY 

C SOLVING AERODYNAMIC OPTIMAL PROBLEM 

C 

c 

c DESCRIPTION OF PARAMETERS 

C 

C TRYRAM (I) THE GO ESS ED VALUES OF RAMDAI 

C, . NTP NT P=1 IF CT IS SPECIFIED 

C 

C SIP .NE, 1 IF CP IS SPECIFIED 

C CTP CTP=CT IF NTP=1 

C CTP=CP IF NTP • NE* 1. 

C NTRY NUMBER OF ELEMENTS GIVEN IN TRYRAM (I) 

C ' NOR LENGTH OF ORAMDA (I) 

C .ORAMDA (I) HYDRODYNAMIC ADVANCE COEFFICIENT AT 

C - THE ZEROS OF T (JRAMDA+1 ,Q) 

C NRAM LENGTH OF TRYRAM (I) 

C 

C 

DIMENSION TRYRAM (NRAM) , ORAMDA (NOR) ,X (5) 

COMMON /DATAl/JRAMDA 
COMMON /DATA6/HUB 
COMMON /DATA11/NCIRCU 
COMMON /DATA12/FM,TM,VF,RAMDAP 
COMMON /DATA14/G (15) 

COMMON /DATA15/CT,CP,ETAI 
COMMON /DATA24/SQHEAN 
COMMON /DAT A3 6/1 ME AN, I PR BLH 
COMMON /DATA37/ENS (2) 

COMMON /DATA38/XX (2,15) 

COMMON /DATA39/GT (2) ,GP{2) 

JX=JRAMDA+1 

IXX=NTRY+1 

40 FORMAT (//,5X, '•************************************««) 
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41 FORMAT (5X , ******* *****«) 

42 FORMAT (5X , ******* OUTPUT FROM AERODYNAMIC *****♦») 

43 FORMAT (5X,»***** OPTIMUM PROBLEM - *****»») 

44 FORMAT (5X, «=M:*************************:********* n) 
WRITE (6,40) 

WRITE (6,41) 

WRITE (6,42) 

WRITE (6 ,43) 

WRITE (6, 41) 

WRITE (6,44) 

DO 51 11=1, IXX 
C 

c' COMPUTE TEE FIRST GUESS OF RAMDA (OQZERO (I) ) 

C THAT IS, ORAMDA(I) 

C 

DO 52 J=1,JX 
ORAMDA (J)=TSYRAM (II) 

52 CONTINUE' 

WRITE (6,56) ORAMDA (1) 

56 FORMAT (//,5X, »RAMDAI=» , F15. 10//) 

CALL RAMCOS (ORAMDA, NOR) 

CALL HATIJ 
CALL CIRCUL 
CALL UATCHBY 

WRITE (6,5) (I,G (I) ,1=1 , MCIRCU) 

5 FORMAT (5X, "G («,I2, M ) =» ,F25, 10) 

CALL AEROCF 

WRITE (6, 6) CT,CP, ETAI 

6 FORMAT (//,5X,"CT=» ,F10«, 5//5X , «CP=” , FIG. 5//,5X, 

S "IDEAL EFFICIENCY ETAI=»,F10, 5//) 

WRITE (6,57) 

57 FORMAT (//,10X,"++R++ ,, ,12X,*'++0++' , /11X,» ++UA++», 

$ 1 1X, " + + UT++» , ,9X, , '++RAMDAI++'», 8X , "++GAMMA++"//) 

R=HUB 

7 IF £R* GT» 1, }' GO TO 9 

Q= (2 a *R- (1« +HUB) )/ (1 *“ HUE) 
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CALL VELOCIT(UA,UT,Q) 

RAH=R* (VF+UA) / (R/RAMDAP+UT) 

WRITE (6,8) R,Q, UA,OT / RAH , CIRLF (Q) 

R=R+G'* 1 
GO TO 7 

8 FORMAT (6 (2X, FI 5, 7) ) 

9 IF (NTEY* EQ„ 1) GO TO 55 
IF (II* EQ* IXX) GO TO 55 
IF(NTP.EQ.I) X(II)=CT 
IF (NTP* NE* 1 ) X{II)=CP 
IF {II.LT. NTRY) GO TO 51 

C 

C INTERPOLATE THE OPTIMUM HYDRODYNAMIC 

C ADVANCE COEFFICIENT 

Y=CTP 
SUM=0. - 

DO 53 1=1 # NTRY 

A=1i : 

B=1, . 

DO 54 J=1,NTRY 
IF (J« EQ* I) GO TO 54 
A=A* (Y-X (J) ) 

3=B* (X (I) -I (J)) 

54 CONTINUE 

SUM=SUH+ (A/B) *TRYRAM (I) 

53 CONTINUE 

TRYRAM (NTRY+1) =SUM 
51 CONTINUE 

55 IF(IMEAN«NE, 1) RETURN 
CALL MEANSQ 

ENS (1)=SQHEAN 
WRITE (6,6C) ENS (1) 

60 FORMAT (//5X # "ENSEMBLE MEAN SQUARE OF»/#5X, 
$ 11 THE ACOUSTIC PRESSURE; »/,5X,E25, 12//) 

IF (IPRBLH • NE. • 4) .RETURN 
- GT{1) =CT 
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GP (1) =CP 

DO 61 J=1,JX 

XX (1,J) =ORA MDA (U) 

61 CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE NQPTIML {OSASDA,NOS,NTP,CTP) 

C 

C PURPOSE 

c 

c TO OBTAIN THE SECOND FEASIBLE SOLUTION OF 

C HYDRODYNAMIC ADVANCE COEFFICIENT BY SOLVING THE 

C NON-OPTIMUM PROBLEM WITH GIVEN CHARACTERISING 

C FUNCTION 

C 

c 

REAL KTP 

DIMENSION ORAMDA (NOS) ,GG (15} ,RAM (15) 

COMMON /DATA1 /JRAMDA 
COMMON /DATA6/HUB 

COMMON /DATA8/NBLADE,KWBNH, LUPIND 
COMMON /CON ST/PAI 
COMMON /DATA11/MCIRCU 
COMMON /DATA1 2/FM, TM , VF f RAMDAP 
COMMON /DATA14/G (15) 

COMMON /DATA15/CT,CP,ETAI 
COMMON /DATAl 6/OQZERQ (15) 

COMMON /DATA35/CT1,CT2,CP1,CP2 
COMMON /BATA24/SQMEAN 
COMMON /DATA36/IMEAN ,IPRBLM 
COMMON /DATA37/ENS (2) 

COMMON /DATA38/XX (2 # 15) 

COMMON /DATA39/GT (2) ,GP (2) 

20 FORMAT (//,5X # »******* ! <'**::fc=fe**#*******: : *****#tt) 
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21 FORMAT (5X, "***** *****»') 

22 FORMAT (5 X, "***** OUTPUT FROM NON- *****■»} 

23 FOBMAT (5X, H ***** OPTIMUM PROBLEM *****«•) 

24 FORMAT (5 X, •’ ****************************” ' ^ 

WRITE {6 ,20) 

WRITE (6,21) 

WRITE (6,22) 

WRITE (6,23) 

WRITE (6, 21) 

WRITE (6,24) 

C 

C*** INPUT.. DATA 
C 

c 

C INPUT THE MAXIMUN NUMBER OF ITERATIONS 

C 

IMAX=3C 

C 

C INPUT THE COEFFICIENTS OF CHARACTERISING FUNCTION 

C 

GG (1) =0*032 67 - 

GG (2) =-0,00083 
GG (3) =-0.00467 
GG (4) =-0. G0243 
GG (5) =-0.00027 
IF (MCIRCU*LE* 5) GO TO 1 
DO 2 M=6 ,MC IRCU 
GG (M) =0. ' 

2 CONTINUE 

1 DO 3 M=1,MCIRCU 
G (M) =GG (M) 

3 CONTINUE 
C 

C 

C*** END OF INPUT' 

C 
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C 

C FIRST APPROXIMATION 

C 

JMAX= JRAMDA+1 
DO 4 J=1,JMAX 
ORAMDA (J) =RAMDAP*VF 

4 CONTINUE 
1=1 

WRITE (6,5)1 

5 FORMAT ( /5 X , H ITER AT 10 N” ,2X,I3/) 

6 CALL RAMCOE (ORAMDA, NOR) 

CALL HATIJ 

CALL UATCHBY 
CALL CTP12 

IF (NTP. HE. 1) GO TO 7 

KTP= (-CT1+SQRT (CTl *CTl +4* *CTP*CT2) )/ 

$ (2**CT2) 

GO TO 8 

7 KTP= (-CP1+SQRT (CP1*CP1 +4.*CTP*CP2)') / 

$ (2* *CP2) 

C 

C COMPUTE NEW RAHDAI AT OQZERO (I) 

C 

ICHK=2 

8 DO 9- J=1 , UMAX 
Q=OQZERO (J) 

R=0,5* ((1,-HUB) *Q+1« +HUB) 

CALL YELOCIT (UA,UT,Q) 

RAM (J) =R* (VF+KTP*UA) /(R/RAMDAP+- 
$ KTP*UT) 

9 CONTINUE 

DO 11 J=1 , JMAX 

IF (ABS (RAM (J)-ORAMDA (J) ) . GT, 1,E-6) GO TO 12 
11 CONTINUE 

DO 13 J=1 , JMAX 

ORAMDA (J) =0*5* (ORAMDA ( J) +RAM { J) ) 
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13 CONTINUE 
GO TO 10 

12 DO 14 J=1,JMAX 

ORAMDA (J) =C.5* (ORAMDA (J) + RAM (J) ) 

14 CONTINUE 
I-I+1 

IF (I. GT* IHAX) GO TO 10 
HR ITS (6,5)1 
WRITE (6,15) KTP 

15 FORMAT (5X,"KTP=",F15. 8/) 

WRITE (6,16) (ORAMDA (J) ,J=1,JMAX) 

16 FORMAT (5X,« VALUE 5 OF RAM DAI" ,/5 (FI 5 , 8 , 2X) ) 

GO TO 6 

10 DO 17 M=1,MCIRCU 
G (H) =KTP*GG (M) 

17 CONTINUE 
CALL BATCHBY 
CT=KTP*CT1+KTP*KTB*CT2 
CP«KTP*CP1+KTP*KTP*CP2 
BTAI=CT*V F/CP 

WRITE (6,35) (I, G (I) ,I=1,SCIRCU) 

35 FORMAT (5X,"G("„I2,») ='\F25, 10) 

WRITE (6,36) CT, CP, ETAI 

36 FORMAT {//,5X,«CT=», FIS. 5//5X,"CP=«, FIS, 5//,5X, 

$ "IDEAL EFFICIENCY ETAI=" ,F1G. 5//) 

WRITS (6,57)- 

57 FORMAT (//,1CX,"++H++», 12X,»++Q++" ,1 1 X , "++UA++", 
$ 1 1 X , "++ GT++ " , 9X , « ++ R A MD AI ++ ", 8X , "++GAMMA++"//) 
R*HDE 

37 IF (R, GT* 1*) GO TO 39 
Q=(2. *R- (1.+HUB) )/ (1,-HUB) 

CALL YELOCIT (UA, UT ,Q) 

RAMD=H* (VF+UA) / (R/RA SDAP+DT) 

WRITE (6,38) R,Q,UA,OT,RAMD,CIRLF (Q) 

R=R+0, 1 

GO TO 37 
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38 FORMAT (6 (2X,F1 5» 7) ) 

3 9 IF (IMEAN* NE« 1) RETURN 
CALL MEANSQ 
ENS (2) =SQMEAN 
WRITE (6,40) ENS (2) 

40 FORMAT (//,5X, “ENSEMBLE MEAN SQUARE OF»/*5X, 
$ “THE ACOUSTIC PRESSURE:”, 2X,E 25* 12//) -- 

IF (IPRBLH , NE. *4) RETURN 
DO 41 J=1,JMAX 
XX (2, J) =ORAMDA (J) 

41 CONTINUE 
GT (2)=CT 
GP (2) =CP 
RETURN 
END 

C 

C 

SUBROUTINE CTPl 2 
C 

C PURPOSE 

C 

C TO COMPUTE CT1, CT2, CPI , AND CP 2 

C 

COMMON /CON ST/PAI 
COMMON /DATA6/HUB 

COMMON /DAT A8/N BLADE , KWENH, LU PIND 
COMMON /DATA11 /MCIRCU 
COMMON /DAT Al 2/FM, TM , VF , RAMDA P 
COMMON /DATA14/G (15) 

COMMON /DATA1 5/CT,CP , ETAI 
COMMON /DATA31 /IVE1DG , IVELCHB 
COMMON /DATA33/UACEB (51) ,UTCHB(51) 

COMMON /DATA35/CT1 ,CT2,CP1,CP2 
NB=NBL ADE 
C 

C CALCULATE CTl AND CPI 
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C 

FACT=G.5*PAI*(1.-HUB) *FLCAT(NB) 

CT1=FACT* {G (1) * (1.+HUB) +G (2) *0,5 * {1,-HOB) )/ 

$ RAMDAP 

CPl=FACT* (G (1) * (1, +HOB) *VF+G (2)* 0,5 *(1,-HOB)* 
$ TF) /RAMDAP 
C 

C COMPUTE CT2 AND CP2 

CT2=0, 

CP2=0 » 

FX=G,5* (lo+HTJB) 

FY-0, 25* (1 b -HUB) 

DO 1 H=1,MCIRCU 

CT2=CT2+G (M ) * {UTCHB (M) -UTCHB (M+2) ) 

CP2=CP2+G (M) * ( (UACHB (M) -OACHB (M+2) ) *FX+ 

$ {UACHB {IABS {M-2) +1) -UACHB {M+3) } *FY) 

1 CONTINUE 

CT2=FACT*CT2 

CP2=FACT*CP2/RAMDAP 

RETURN 

END 

D, 

C 

SUBROUTINE AEROCF 
C 
C 

C PURPOSE 

C 

C TO COMPUTE THE THRUST AND P09ER COEFFICIENTS 

C 

C 

COMMON /CCNST/PAI 
COMMON /DATA6/HUB 

COMMON /DATA8/NBLADE, KWENH, LUPIND 

COMMON /DATA11/MCIRCU 

COMMON /DATA12/FM r TM, VF, EAMDAP 
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COMMON /DATA14/G (15) 

COMMON /DATA15/CT,CP,ETAI 
COMMON /DAT A31/IVELDG, IVELCHB 
COMMON /DATA33/OACHB (51) r DTCHB (51) 
NB=NBLADE 

FACT-C *5*PAI* (1. -HUB) AFLOAT (NB) 

FX=G* 5*(1. +HOB) 

FY=0.25* (1.-HUB) 

SDMT=G (1)* (1 » +HU3) /RAM DAP 
SUMT=SUMT+G (2)* (1, -HUB ) / (2. *BAMDAP) 

SUMP=G (1) * (1. +HUB) *VF 
SUMP=SUMP+G (2) * (1* -HUB) *VF/2. 

DO 1 M=1 , HCIRCU 

SUHT=SUMT+G (H) * (UTCHB (M) -UTCHB (M + 2) ) 
SUMP=SUMP+G (M) * (FX* (UACHB (M) — UACHB (H+2) ) + 
$ FY* (UACHB (IABS (M-2) +1) -UACHB (M+3) ) ) 

1 CONTINUE 

CT=SUHT*FACT • 

CP=SUMP*FACT/RAMDAP 

ETAI=CT*VF/CP 

RETURN 

END 

C 

C 

FUNCTION RAIJKM(I,J,K,M) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE RAMAD (I) (J) (K) (M) 

C 

C 

COMMON /CON ST/PAI 

ISUM=0 

MJ=M-J 

KI=K-I 
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IABSMJ=IABS (MJ) 

IABSKI=IABS (KI) 

IF (MJ« HE* 0) M J=M J/IABSMJ 
IF(KI.NE.O) KI=KI/IABSKI . 

IF ( (I +K-J-M) • EQ« 0) ISUM=1 
IF { (I+K-IABSMJ) * EQo 0) ISUa=ISUH+MJ 
IF ( (IABSKI-J-H) • EQ® 0) ISUM=ISUM+KI 
IF (IABSKI. EQ. IABSMJ} ISUM=ISUM+KI*MJ 
BAIJKM=PAI* FLOAT (ISOM) /8, 

RETURN 

END 

C 

c 

SUBROUTINE CIRCUL 
C 
C 

C PURPOSE 

C 

C TO COMPUTE THE COEFFICIENTS OF THE CHEBYSHE? 

C EXPANSION OF THE CIRCULATION 

C 

C CIRCULATION (Q) =SQRT ( 1-Q**2) *SUMMATION G(M)* 

c U(a-1,Q), a=l(l)aciRcu 

C 
C 

C DESCRIPTION OF PARAHETEKS 

C 

C NT LENGTH OF TF 

C, NU LENGTH OF UF 

C 
C 

• COMMON /CONST/PAI 

COMMON /DATA3/IDEG, ICHEBY 
COMMON /DATA4/HA (31,31) ,HT (31,31) 

COMMON /DATA6/HUB 
COMMON /DATA1 0/QCIRCU (15) 



142 


COMMON /DATA 11 /MCIRCU 
COMMON /DATA12/FM,TH,VF,RAMDAP 
COMMON /DATA14/G (15) 

DIMENSION A (15,15) ,B (15,1) , RAMDAL (15) 
DIMENSION WK (300) ,TF (31),aF(31) 

C COMPUTE RAMDA AT QCIRCU(L) 

DO 1 1=1 , MCIRCU 
Q=QCIRCU (L) 

RAM DAI (L)=RAKDAF (Q) 

1 CONTINUE 

C CONSTRUCT TEE COLUMN VECTOR B 

DO 2- 1=1, MCIRCU 
B (L) = RAMDAL (L) /RAM DAP- VF 

2 CONTINUE 

C ' CONSTRUCT THE COEFFICIENT MATRIX A 

IF (MCIRCU. GE. IDEG) IMAX=MCISCU 
IF (MCIRCU. IE. ID EG) IMAX=IDEG 
DO 3 L=1 , MCIRCU 
Q=QCIRCU (I) 

C 

C INPUT DATA 

C 

NT=31 

NU=31 

C 

CALL TCHEBY (TF , IMAX, Q, NT) 

CALL UCHEBY (UF, IMAX, Q, NT) 

FA=2.*RAHDA1 (L) /( (1. -HUB)*QCIRCU (L) + 1.+HUB) 
DO 5 H=1, MCIRCU 
FAAM=PAI*FLOAT (M) / (1 • -HUB) 

SUM=0 3 

1=0 

6 HI=1 

IF (I, GT» IDEG) GO TO 7 
IF(I.SQ.O) HI=2. 

SA=0. ' 
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J=G 

8 HJ=1 » 

IP (J. ST, IDE G) GO TO 9 
IF(J, EQ, 0) HJ=2. 

IF(J. GT..H) GO TO 10 

SA=SA+ (HI (I+1,J + 1) -FA*HT (1+1 ,J+1) ) *TF (1 + 1) * 

$ TF (J+1)*UF (H) *FAAM/ (HI*HJ.) 

J=J+1 
GO TO 8 

10 SA=SA+ (HA (1+1, J+1) -FA*HT (1 + 1 , J+ 1 ) ) *TF (1+1) * 

$ TF(M + 1)*UF (J) +FAAM/ (HI*HJ) 

J=J+1 
GO TO 8 

9 SUM=SUM+SA 
1=1 + 1 

GO TO 6 
7 A (L,M) =SIJH 
'5 CONTINUE 
3 CONTINUE 
C 

C SOLVE THE SYSTEM OF EQUATIONS A (L, M) *G (M) =B (L) 

C FOR G (M) ,L,M=1 (1) KCIBCU 

C BY USING THE LIBRARY SUBROUTINE LEQT2F OF U OF I 

C 

C 

IER-500 

IDGT=0 

H=1 

N=HCISCU 

IA=15 

CALL LEQT2F (A,a,N,IA,B,IDGT,WK,IER) 

DO 20 H=1,aCIRCU 
G (H)*B (H, 1) 

20 CONTINUE 

IF { { (IEH-34) * (IER-12 9) * (IER-131) ) , NE.,0) RETURN 
WRITE (6,24). 
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24 FORMAT (///, 5X, ***** ITERATION IS TERMINATED "/t 
$ 5X,»DUE TO ILL-CONDITIONED MATRIX ***") 

STOP 

END 

C 

C 

SUBROUTINE UATCHBY 
C 

C PURPOSE 

C 

C TO COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV 

C EXPANSIONS OF THE INDUCED VELOCITY COMPONENTS 

C 
C 

COMMON /DATA31 /IVELDG, I VELCHB 
COMMON /DAT A32/QVEL ( 51 ) 

COMMON /DAIA33/UACHB (51) ,UTCHB (51) 

DIMENSION A (50) ,.UA (51), UT (51) 

C 

C COMPUTE THE INDUCED VELOCITY COMPONENTS 

C 

JMAX=IVELCHB+1 
DO 1 J=1 , JMAX 
Q=QVEL (J) 

CALL INDVEL (Q, AN SA , ANST) 

UA (J) =ANSA 
UT(J)=ANST 
1 CONTINUE 
C 

C COMPUTE UACHB (I) ,AND UTCHB (I) ,1=1 (1 ) (IVELDG+1) 

C 

C NA=LENGTH OF ARRAY A 

C 

NA=5G 

NLETH=51 

MCHEBY=JMAX 
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NDEG^IVELDG 

CALL CHEBCF (UA , HCHEB Y, NDEG ,NLETH # AO, A, NA) 

UACHB (1) =A0 
DO 2 J=1 , NDEG 
UACHB (J+1 ) “A (J) 

2 CONTINUE 

CALL CHEBCF {UT,MCH2BY, HDEG , NLETH, AO, A, BA) 

DO 3 J=1,NDEG 
UTCHB (J+1) -A (J) 

3 CONTINUE 
UTCHB (1) -AG 
RETURN 

END 

e 

c 

<■* 

C 

SUBROUTINE HATIJ 
C 
C 

C PURPOSE 

C 

C TO COMPUTE THE COEFFICIENTS OF THE DOUBLE CH2BYSHE? 

C EXPANSION OF THE AXIAL AND TANGENTIAL INDUCTION 

C FACTORS, HA (-1 , J) AND HT(I,J) ,THAT IS, 

C 

C IA {Q,QO)=SUHHATION** HA (I, J) *T (I,Q) *T(J,QO), 

C I, J-3 (1) IDEG 

C 

C IT (Q,QO) =SUMj3ATION * 1 HT (I, J) *T (I ,Q) *T(J,QO), 

C I,J— 0 (1) IDEG 

C 

C DESCRIPTION OF PARAMETERS 

C 

C QIND(I) THE ITH ZERO OF T (ICHEBI+1 ,Q} 

C 



QIND(I) =COS ( (2*1+1 ) *P AI/ (2* (ICHEBY+1) ) 
1=0 (1) ICHEEY 


COMMON /DATA3/IDEG,ICHEBY 
COMMON /DAT A 4/HA (31 ,31) , BT (31,31) 

COMMON /DATA5/QIND (51) 

COMMON /DATA6/HUB 

.COMMON /DATA8/NBLA DE,KWENH,LUPIND 

COMMON /DATA9/EP SIND 

DIMENSION AXIAL (51 ,51) , TANG (51,51) 

DIMENSION S (51) , BO (51) ,R AM DA (51) 

DIMENSION PRJ (51 ,51) , HAT (51,51) 

COORDINATE THAN SFOEM ATION R= ( (1-HUB) *Q+1 + HDB)/2 
IX-1 +ICHEBY 
DO 1 1=1, IX 

R (I) =0« 5 * ( (1,-HUB) *QIND (I) +1*+HUB) 

HO (I) =R (I) 

CONTINUE 

COMPUTE RAM DA AT E 
DO 2. 1=1, IX 
Q=QIND (I) 

EAMDA (I)=BAMDAF(Q) 

CONTINUE 

COMPUTE THE AXIAL AND TANGENTIAL INDUCTION 

FACTORS 

NB=NBLADE 

NPR=IX 

NPRO=IX 

INPUT DATA 

KR, KHO : DIMENSION OF AXIAL, TANG 

KR=51 

KRO=51 

KCONT=KWENH 

EPSLON=EPSIND 
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LUPER=LUPIND 

CALL INDFACT (NB, NPR, NPRO ,R, RO , RAMDA , AXIAL, . 

$ TANG,KR,KRO,KCONT,EPSLON,LUPER) 

C COMPUTE HA (I, J), AND HT (I,J) 

NC=KR 

•x 

C 

C NA: THE LENGTH OF HA (NA, NA) AND HT(NA,NA) 

NCHEBY=ICHEBY 
NDEG=IDEG 
NAT-51 
NHAT=IDEG+1 

CALL DOOCHB (AXIAL, NCHEBY,ND3G, HAT, NAT, NC, PRJ, R) 
DO 3 1=1, NHAT 
DO 4 J=1 , NHAT 
HA (I,J)=HAT (I, J) 

4 CONTINUE 
3 CONTINUE 

CALL DOUCHB (TANG,NCHEB Y, NDSG ,HAT , NAT, NC, PRJ, R) 
DO 5 1=1, NHAT 
DO 6 J=1,NHAT 
HT (I, J)=HAT (I, J) 

6 CONTINUE 

5 CONTINUE 
RETURN 
END - 

C 

C 

FUNCTION RAMDAF (Q) 

C 

C 

C PURPOSE 

C 

C TO COHPUTE RAMDA AT Q 

C 

C 

C THAT IS, 



148 


C 

C RAMDAF (Q) =A0RAMD/2 + ARAMDA (J) * T (J,Q) J=1, JRAMDA 

C 

COMMON /DATAl /JRAMDA 

COMMON /DATA2/A OR AMD, ARAMDA (14) 

DIMENSION A (14) 

NLETH=14 
N= JRAMDA 
AG=AGRAMD 
DO 1 I=1,JRAMDA 
A (I) -ARAMDA (I) 

1 CONTINUE 

CALL SUMCHB (AQ ,A,NLETH,N,Q,ANS) 

RAMDAF=ANS 

RETURN 

END 

C 

c 

SUBROUTINE RAMCOE (ORAMDA ,NOS) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV 

C EXPANSION OF SAM DA, .THAT IS 

C 

C RAMDA (Q) =A0RAMD/2 +ARAHD A' ( J) *T (J , Q) J=1 , JRAMDA. 

C 

C DESCRIPTION OF PARAMETERS 

C 

C 

C ORAMDA VALUES OF RAMDA A ZEROS OF 

C T (JHAM DA+1 ,Q) 

C NOR LENGTH OF ORAMDA 

COMMON /CON ST/PAI 
COMMON /DAT Al/JRAMDA 
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COMMON /DATA2/AGRAMD,ARAMDA (14) 

M= JR AM DA 

H=PAI/FLOAT (2* (JRAMDA+ 1) ) 

FACT=2. /FLOAT (JSAMDA+1) 

J=0 

1 IF (J* GT» J1AMDA) GO TO 2 
X=COS (FLOAT (J) *H) 

CALL SUMODD (ORAMDA,M,NORVX,ANS) 

IF(J. EQ.O) AORAMD=FACT*ANS 
IF (J» NE» 0) ARAMDA (J) =FACT*ANS 
J=J+1 
GO TO 1 

2 RETURN 
END 

C 

C 

SUBROUTINE DOUCHB (C, NC HEBI f NDEG, A , NA, NC , PRJ, WK) 

C 

C 

C PURPOSE 

c 

C 

C TO COMPUTE THE COEFFICIENTS OF A. FINITE 

C BIVARIATE CHEBYSHEV EXPANSION 

C THAT IS, 

C 

C F (X,Y) = B (I, J) * T(I,X) * T(J,X), I,J=0 (1) NDEG 

C 

C WHERE T (X,I) =CHEBYSHEV POLYNOMIAL OF DEGREE I 

C 

C B (0,0) =A (1 ,1) /4 

C B (1,0) =A (1+1 ,1)/2 (I • GT» 1) 

C B (0 , J) =A (1 , J+1 ) /2 (J • G T* 1 ) 

C B(I,J)=A(I + 1,J+1) (I,J .GT, 1) 

C 

C 
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C DESCRIPTION OF PARAMETERS 

C . 

C NC 

C NA 

c c 

c 
c 
c 

n 

Var 

c 

C PRJ 

C WK 

C 

DIMENSION C (NC , NC) ,A (NA,NA) , PRJ ( NC , NA) ,WK(NC) 
COMMON /CONST/PAI 
H=PAI/FLOAT (2* (NCHEBY+1) } 

IR MAX— NCHEBY+1 
IS MAX= NCHEBY+1 
IMAX=NDEG+1 
JMAX=NDEG+1 
C COMPUTE PRJ 

DO 1 JX=1,JMAX 
J=JX- 1 

X=COS ( H* FLOAT (J) } 

DO 2 IRX=1,.IRHAX 
C CONSTRUCT WK 

DO 3 ISX=1,ISMAX 
WK (ISX)=C (IRX ,ISX) 

3 CONTINUE 
M=NCHEBY 

CALL SUMODD (WK,H,NC,X, ANS) 

PRJ (IRX,JX)=ANS 
2 CONTINUE 
1 CONTINUE 

C COMPUTE COEFFICIENTS A(I,J) 

FACT= A*/ ( (NCHEBY+1) **2) 


LENGTH OF C 
LENGTH OF A 

FUNCTIONAL VALUES OF F { X/ Y) AT THE 
DESCRETE SET (NCHEBY + 1) (NCHEBY+1) POINTS 
(X(I),Y(J)), I, J=Q (1) NCHEBY 
' X (I) =Y (I) THE ROOTS OF T (NCHEBY+1 ,X) 

X (I) =COS ( (2*1) *PAI/ [28 (NCHEBY+1)’) , 

I=G (1) NCHEBY 

WORK AREA (DIMENSION (NC,NA) 

WORK AREA (DIMENSION (NC) ). 
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DO 4 11=1, IBAX 
I=IX“1 

X=COS (H*FLOAT (I) J 
DO 5 JX=1,JMAX 
C. CONSTRUCT SK 

DO 6 IRX=1,IRHAX 
WK (IRX) = PR J (IRX , JX ) 

6 CONTINUE 
- M=NCHEBY 

CALL SUMODD (WK, M, NC, X , ANS) 

A {IX, JX)=FACT*ANS 
5 CONTINUE 
4 CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE CHEBCP { FTBI ,MCHEBY ,ND3G , NLETH , AO , A, NA) 

C 

C 

C - PURPOSE 

C 

C 

C TO EVALUATE " THE' COEFFICIENTS OF THE CHEBYSHEV 

C EXPANSION OF THE FORM F(X)=A0/2 + A (K) * T(K,X) 

C. USING CHEBYSHEV QUADRATURE 

C 

C DESCRIPTION OF PARAMETERS 

C 

C MCHEBY THE NUMBER OF POINTS USED IN 

C CHEBYSHEV QUADRATURE FORMULA 

C FTHI FUNCTIONAL .VALUES AT THE SOOTS OF 

C T (MCHEBY, X) ,THI (I)=COS ( (2*1-1) *PAI/ 

C (2*MCHEBY) ) , 1=1 (1 ) MCHEBY 

C NDEG DEGREE OF CHEBYSHEV EXPANSION 

C NLETH LENTH OF THI, ETA, FTHI, FETA 
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C NA LENGTH CF A 

C AO, A COEFFICIENTS OF CHEBYSHEV EXPANSION 

C 

C 

DIMENSION A (NA) ,FTHI (NLETH) 

COMMON /CONST/PAI 
HODD=HCHEBY-1 
H=PAI/FLOAT (2*HCEBBY) 

K=0 

1 IF (K« GT* NDEG) GO TO 2 

C COMPOTE ALPHAK 

X=C0S (H*FLO AT (K) ) 

CALL SUMODD (FTHI ,MODD, NLETH, X, ANS) 

ALPH AK=2* *ANS /FLOAT (MC HEBY) 

IF (K. EQ. 0) A0=ALPH AK 
IF(K.NE.O) A (K) =ALPHAK 
K=K+1 
GO TO 1 

2 RETURN 
END 

C 

C 

SUBROUTINE SUM0DD{C, H, NLETH, I, ANS) 

C 

c 

C PURPOSE 

C 

C TO COMPUTE THE ODD POLYNOMIALS IN 

C CHEBYSHEV FORM iT X 

C THAT IS, 

C 

C ANS=P(X)=C(J+1)*T{2*J+1,X) , J=Q(1)H 

C 

C DESCRIPTION OF PARAMETERS 

C 

C NLETH LENGTH OF C 
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C 

DIMENSION C (NLETH) 

K=M 

BK1=0„ . 

BK2=0a 

I=2.*X*X-1. 

1 BK=2. *T*BKl-BK2+C{K+1) 

IF (K*EQ,G) GO TO 3 

K=K-1 

BK2=BK1 

BK1=BK 

GO TO 1 

3 ANS=X* (BK-BK1) 

RETURN 

END 

C 

C 

SUBROUTINE SUMCHB (AO , A , NLETH, N ,X,ANS) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE VALUE OF POLYNOMIALS IN 

C CHEBYSHEV FOBHS 

C THAT IS, 

C 

C ANS=P (X)=A0/2 + A (K) *T (K,X) , K=1 (1 ) N 

- > 

C 

C DESCRIPTION OF PARAMETERS 

C 

C NLETH LENGTH 0 A 

C 

DIMENSION A (NLETH) 

BK1=G, 

BK2=0, 

K=N 
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1 IF (K« EQ» 0) AK=AO 
IF (K* HE a 0) AK=A(K) 

BK=2» *X*BKl-BK2+AK 
IF (K« EQo 0) GO TO 2 
K=K-1 

BK2=BKl 
BKl=BK 
GO TO 1 

2 ANS=0. 5* (BK-BK2) 

RETURN 

END 

C 

c 

SUBROUTINE ZEROS (NDEGRE, NZERO, ZE RO) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE ZEROS OF T {NDEGRE, Q) 

C 

C DESCRIPTION OF PARAMETERS 

C 

C ZERO (I) ZERO (I) =COS ( (2*1-1) *PAI/ (2*NDEGRE) ) 

C 1=1(1) NDEGRE 

C 

C NZERO LENGTH OF ZERO 

C 

C 

COMMON /CONST/PAI 
DIMENSION ZERO (NZERO) 

H=PAI/F10AT (2*N DEGRE) 

DO 1 1=1, NDEGRE 
Z=H*FLOAT (2*1-1) 

ZERO (I) =COS (Z) 

1 CONTINUE 
RETURN 



END 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE DOUSUH (A, B, M, N, NM , NN,X , I , ANS) 


PURPOSE 

TO COMPUTE THE SUM OP A FINITE DOUBLE 
CHEBYSHEV S3RIES*,THAT IS, 

ANS=SUMMATIGN* » B (I, J) * T(I,X) * T(J,Y) 

i=G {1 ) a, j=0(1)n 

WHERE B (I,J)=A (1+1, J+1) 

DESCRIPTION OF PARAMETERS 

A COEFFICIENTS OF THE DOUBLE 

CHEBYSHEV SERIES 
B WORK AREA OF LENGTH NM 

N DEGREE IN X 

H DEGREE IN Y 

.NM LENGTH OF THE ROW OF A 

NN LENGTH OF THE COLUMN OF A 


DIMENSION A (NM,NN) ,B (NM) 
COMPUTE B (I) 

IMAX=M+1 

JMAX=N+1 

DO 1 IX=1,IMAX 

I=IX-1 

J=N 
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DNl=0. 

DN2=0. 

2 DN=A(I+1, J+1) +2. *Y*DNl-DN2 
IF(J.EQ.O) GO TO 3 

J=J-1 
BN2=DN1 
DNl=DN 
GO TO 2 

3 B (1+1) =0.5 * (DN-DN2) 

1 CONTINUE 

I=M 
G 111 "=0 ■» 

GM2=0» 

4 GM=B{I+1)+2,*X*GMl-GM2 
IF (I. EQ. 0) GO TO 5 
1 = 1-1 

GM2=GMl 
GHl=GM ' 

GO TO 4 

5 ANS=0 a 5 * (GM-GM2) 

RETURN 

END 

C 

C 

SUBROUTINE FACTOR (Q, QO , AXIAL, TANG) 

n 

C PURPOSE 

c 

C TO COMPUTE THE AXIAL AND TANGENTIAL 

C INDUCTION FACTORS BY USING THE 

C DOUBLE CHEBYSHEV EXPANSIONS 

C 
C 

COMMON /DATA3/IDSG,ICHEBY 
COMMON /DAT A 4/HA (31,31) , HT (31,31) 
DIMENSION A (31,31) ,B (31) 
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NM-31 

NN=31 

H=IDEG 

N=IDEG 

IMAX=IDEG+1 

DO 1 1=1 , IM AX 

DO 2 J=1,IMAX 

A(i r j)=aa (i # j) 

2 CONTINUE 
1 CONTINUE 

CALL DOUSUa (A,B,a,N,.Na,NN,Q,QO,ANS) 

AXIAL=ANS 
DO -3 1=1, MAX 
DO 4 J=1,IHAX 
A (I, J) =HT (I,J) 

4 CONTINUE 

3 CONTINUE 

CALL DOUSUa (A,B,H,N, NH # NN, Q, QO, ANS) 

TANG=ANS 
EE TURN 
END 
C 
C 

FUNCTION CIHLF(Q) 

C 

C 

C PURPOSE 

C 

C TO CONFUTE THE CIRCULATION AT Q 

C ' 

C CIRLE (Q) =SQRT (1-Q**2) * SUMMATION G {M) * U(M-1,Q) 

C 

C 

COHNON /DATA11/HCIRCU 
COMMON /DATA14/G (15) 

N=HCIRCU-1 
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BKl=Oe . 

BK2=0» ' 

K=N 

1 BK=G (K+1 ) +2, *Q*BKl-BK2 
IF (K, EQ.1) GO TO 2 
K=K-1 

BK2=BK1 
BK1=BK 
GO TO 1 

2 ANS^G (1) -BK1+2* *Q*BK 
CIRLF=SQRT(1„-Q*Q) * ANS 
RETURN 

END 

C 

C 

SUBROUTINE INDVEL (Q, UA, UT) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE AXIAL AND TANGENTIAL 

C COMPON ENTS OF THE INDUCED VELOCITY 

C 
C 

COMMON /CONST/PAI 
COMMON /DATA3/IDEG,ICHEBY 

COMMON /DATA4/HA {3 1 , 31 ) , HT (31 , 31 ) 

■> 

COMMON /DATA6/HUB 
COMMON /DATA11/MCIRCU 
COMMON /DAT A1 4/G {1 5) 

DIMENSION TF (31) ,tJF(31) 

I-F (MCIRCU* GE» IDEG) IHAZ=MCIRCU 
IF (MCIRCU, LE, IDEG} IMAX=IDEG 
NT=31 
NU=31 

CALL TCHEBY (TF, IMAX , Q, NT) 
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CALL UCHEBY (UF,IMAX, Q, NU) 

SUMT=0. 

SUHA=0« 

DO 2 M=1,HCIBCU 

FH=G (E) *PAI*FLOAT ( 33} / (I. -HUB) 

1=0 

S UMAI=0<» - 

SUMTI=0. 

\ 

3 HI=1 « ' 

IF (I, GT.IDEG) GO TO 4 
IF (I* EQ* 0) HI=2. . 

J=0 

SUHAJ=G» 

SUMTJ=0, 

5. HJ=1 • 

IF(J, GT.IDEG) GO TO 6 
IF (J. EQ* 0) HJ=2. 

IF ( J. GT. S) GO TO 7 

SUSAJ=SUHAJ+HA [1+1 ,J+1) *TF (1 + 1) *TF (J + 1) * 
$ UF (H) / (HI*HJ) 

SUHTJ=SUaiJ+HT (1+1 ,J + 1) *TF (1 + 1) *TF (J+1) * 
$ UF (M) / (HI* H J) 

J=J+1 
GO TO 5 

7 SUHAJ=SUHAJ+HA (1+1 , J+ 1 ) *TF (1+1) *TF (3+1) * 
$ UF (J) /(HI*HJ) 

SUMTJ=SUHTJ+HT (1+1 , J+1 ) *TF (1 + 1) *TF (H+1) * 
$ UF (J) /(HI*HJ) 

J=J+1 
GO TO 5 

6 SUMAI=SUaAI+SUEAJ 

susti=suhti +suhtj 
1=1+1 ' 

GO TO 3 

4 SUHA=SUHA+SUMAI*FH 
SUHT=SUHT+SUHTI*FH 
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2 CONTINUE 
UA=SUMA 
BT=SUMT 
RETURN 
END 
C 
C 

SUBROUTINE YELOCIT (U A,UT ,Q) 

C 

C PURPOSE 

C 

C TO INTERPOLATE THE INDUCED VELOCITY COMPONENTS 

C BY EVALUATING THE CHEBYSHEV EXPANSIONS OF 

C THE INDUCED VELOCITY COMPONENTS 

C 

COMMON /DAT A31 /I VELD G, IV ELCHB 
COMMON /DATA33/OACHB (51) ,UTCHB (51) ' 

DIMENSION A (50) 

N1ETH=50 
N=IVELDG 
DO 1 1=1 ,N 
A (I)=UACHB (1+1) 

1 CONTINUE 
A0=UACHB (1) 

CALL SUMCHB{AO,A,NLETH,N,Q,ANS) 

UA=ANS 
DO 2 1=1, N 
A (I)=UTCHB (1+1) 

2 CONTINUE 

A0 = UTCHB (1) 

CALL SUMCHB {AO , A ,NLETH,N , Q, ANS) 

UT=ANS 

RETURN 

END 

C 

C 
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C 

SUBROUTINE TCHEBY (T, N, X , NT) 

C 

c 

C PURPOSE 

C 

C TO COMPUTE THE VALUES OF CHEBYSHEV POLYNOMIALS 

C OF THE FIRST KIND OF DEGREES FROM Q-N AT X 

C 
C 

DOUBLE PRECISION Q,TN1 , TN2, TN 
DIMENSION T (NT) 

Q=X 

TH2=1 » DO 
TN1=Q 
T (1 ) =TN2 
T(2)=TN1 

IFfN.LE. 1) RETURN 
IHAX=N+1 
DO 1 I=3,IMAX 
TN=2. DO * Q * TNI - TN2 
T (I) =TN 
TN2=TN1 
TN1=TN 
1 CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE UCHEBY (U,N,X ,NU) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE VALUES OF CHEBYSHEV POLYNOMIALS 

C OF THE SECOND KIND, U(I,X), 1=0 (1)N 
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C 

C 

DOUBLE PRECISION Q,UN,UN1,UN2 
DIMENSION 0 (NU) 

Q=X 

UN2=1 * DO 
UN 1=2, DO * Q 
U (1 ) =UN2 
U(2)=UNl 

IP (N* LE» 1 ) RETURN 
IMAX=N+1 
DO 1 1=3, IM AX 
UN=2,DG * Q* UNI -UN2 
U{I)=UN 
UN2=UNl 
UN1=UN 
1 CONTINUE 
RETURN 
END 
C 
C 

c 

SUBROUTINE INDPACT {NB, NPR, NPRO, H, RO, RAMDA, AXIAL , 
$ TANG,KR, KRO,KCONT,EPSLON,LUPER) . 

C 

C PURPOSE 

C 

C TO COMPUTE TABLES OF THE AXIAL AND TANGENTIAL 

C INDUCTION FACTORS 

C 
C 

C DATA FILE NEEDED 

C 

C COSGRA 

C 

C DESCRIPTION OF PARAMETERS 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NB NUMBER OP PROPELLER BLADES 

NPR NUMBER OF ELEMENTS OF ARRAY R 

NPBO NUMBER OF ELEMENTS OF A EBAY BO 

KR LENGTH OF' ARRAY R 

KRO LENGTH OF ARRAY HO 

AXIAL (K, I) AXIAL INDUCTION FACTOR AT (R(K),RO(I)) 

TANG (K ,1) TANGENTIAL INDUCTION FACTOR AT (R(K),EO(I)) 

1UPER UPPER LIMIT OF IA2 

EPSLON UPPER LIMIT OF IA3, EPSLON « GE« 0 

R ARRAY CONTAINING THE RADIAL COORDINATES 

OF REFERRENCE POINTS 

RO ARRAY CONTAINING THE RADIAL COORDINATES 

OF SOURCE POINTS 

R AND RO_ AREN AS SU MED TO BE EQUAL IF 
ABS {R— RO) ,LE. 1.E-1Q 

RAMDA {K> HYDRODYNAMIC ADVANCE COEFFICIENT AT RO (K) 

KCONT CONTROL PARAMETER 

KCONT « EQ a .1 IF WRENCH MODIFIED FORMULA IS ■ 
USED 

KCONT . NE a . 1 IF ALTERNATE METHOD 

IS USED 


REAL R (KR) ,RO(KRO) , AXIAL (KR, KRO) , TANG (KR, KRO) 
REAL RAMDA (KRO) 

IF (KCONT . NE. -1} GO TO 5 

COMPUTE INDUCTION FACTORS USING WRENCH'S 
MODIFIED FORMULAS 
DO 6 K=1,NPRO 
ROU=RO (K) 

EAMDAK=RAMDA (K) 

DO 7 J=1,NPS 
R3=R(J) 

CALL HRENF (NB, RR, ROU ,RAM DA K, ANSA, ANST) 

AXIAL (J,K) = ANSA 
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TANG (J , K) =ANST 

7 CONTINUE 
6 CONTINUE 

RETURN 

5 DO 1 1=1 , NPR 
DO 2 J=1,NPRO 
AXIAL (I,J}=Q. 

2 CONTINUE 
1 CONTINUE 

CALL INDA12 (EPSLON, NB, LUPER, NPR,NPBO,R,RO,RAMDA, 

$ AXIAL, KR, KRO) 

IF (EPSLQN . LE« 1.E-10) GO TO 8 

CALL INDA3 (EPSLON,NPR, N P RO ,.R , RO , RA MDA , , 

$ AXIAL, KR, KRO) 

8 CALL INDA4 (NPR, NPRO, R,RO,RAHDA, AXIAL, KR, KRO) 

B= FLOAT (NB) 

DO 3 K= 1 , N P RO 
ROU=EO (K) 

RAHDAK=EAMDA (K) 

DO 4 J=1,NPR 
RR=R(J) 

IF (ABS (RR-ROU) .LS. 1.E-10) AXIAL (J,K) =ROU/ 

$ SQRT (RAHDAK*RA 8DAK+RR*R0U) 

TANG (J,K)=RAMDAK* ( (ROU-RR) *B/RAHDAK- 
$ AXIAL (J,K) )/RR 
4 CONTINUE 

3 CONTINUE 
RETURN 
END 

C 

C 

SUBROUTINE INDAl 2 (EPS,NB,LUPER,NPR, NPRO,R, RO,RANDA , 
$ AXIAL, KR, KRO) 

C 
. C 
C 
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C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PURPOSE 

TO COMPUTE 1412 BY USING 4 25-POINT 
GAUSS-LEGENDEE FORMULA ON EACH SUB-SUBINTERVAL 

PRECISION: DOUBLE PRECISION 


DESCRIPTION OF PARAMETERS 


NBLADE 

NLAST 
IB EG IN 

I FINAL 

EPSLON 

L 

ALX 

AUX 

NHN 

NST 

NFI 


NUMBER OF PROPELLER BLADES-DEFINED IN 
MAIN PROGRAM (USES SUPPLIED) 

NUMBER OF SUB^INTESVALS 
SUB-INTERVAL PARAMETER-DEFINED IN MAIN 
PROGRAM (USES SUPPLIED) 

SOB-INTERV.AL PARAMETER— DEFINED IN MAIN 
PROGRAM (USER SUPPLIED) 

E.G. . IBEGIN=1, AND IFINAL=NLAST IF ALL 
ALL SUN-INTERVALS ARE CONSIDERED 
LOWER LIMIT OF INTEGRAL IA12 
2 * PAI * L= OPER LIMIT OF INTEGRAL IA12 
-DEFINED IN MAIN PROGRAM (SUER 
SUPPLIED) 

LOWER LIMIT OF EACH SUB-INTERVAL- 
DEFINED IN MAIN PROGRAM (USER SUPPLIED) 
UPER LIMIT OF EACH SUB-INTERVAL 
-DEFINED IN MAIN PROGRAM (USER SUPPLIED) 
NUMBER OF SUB-SUBINTERVALS OF EACH 
SUB-INTERVAL (ALX (I) , AUX (I) ) 

BLADE PARAMETER-DEFINED IN MAIN PROGRAM 
(USER SUPPLIED) 

BLADE PARAMETER-DEFINED IN MAIN PROGRAM 
(USER SUPPLIED) 
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C 

C 

r 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


E« G* - NST (I) =1 AND NFI(I)=1 IP ALL BLADES ARE 
CONSIDERED IN SOB- INTERVAL- I 
ICHK IF ICHK=1 , THIS PROGRAM CALCULATES IA12 WHICH 

HAS INDEX J» EQ> K AND (I- J) . LE. 1 AND PRINTS OUT 
THESE RESULTS (SEE REMARK) 


REMARK 

THE PARAMETER NHN (I) FOR EACH SUB-INTERVAL MAY BE 
CHOSEN AS FOLLOWS TO OBTAIN DESIRED NUMERICAL ACCURACY 


STEP 1: 
STEP 2: 
STEP 3: 
STEP 4: 
STEP 5: 
STEP 6: 

STEP 7: 


SET ICHK=1 

ASSUME NHN (I) ,I=1,NLAST 
RUN THIS PROGRAM 

INCREASE OR DECREASE NHN (I) ,1=1 , NLAST 
RUN THIS PROGRAM 

COMPARE CURRENT OUTPUT WITH ALL PREVIOUS 
OUTPUT 

REPEAT STEPS 4 TO 6 UNTIL DESIRED ACCURACY 
IS OBTAINED 


IMPLICIT DOUBLE PRECISION (A-H ,0- Z) 

DOUBLE PRECISION ALX (20) , AUX ( 20) 

REAL R (KR) , RO (KRO) , AXIAL (KR, KRO) ,RAMDA (KRO) 
REAL EPS 

INTEGER NST (20) ,NFI(20) ,NHN (20) 

COMMON /BLADE/CB (10) ,SB (10) 

COMMON /ADVAN/RAMDAK 
COMMON /COO SD/SR , RO U 
COMMON /GQUZ5/NH , NST ART , NFINAL 
COMMON /GQUZ6/HLETH 

COMMON /GQUZ7/ZZ (25,5) ,CZZ- (25 ,5) , SZZ (.25, 5) 
EXTERNAL FAXIAL 

PAI=3o 1 41 59 2653 5 8979323 8 46 26 433832 7 9D0 
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NBLADE=NB 

C 

C*** INPUT DATA 

C 

ICHK=2 

NLAST=6 

EPSLON=EPS 

IBEGIN=1 

IFINA1=NLAST 

L=LUP3E 

ALX (1) =EPSLON 

AUX (1 } =PAI/1 00 » DO 

NHN (1) =1 

NST (1 ) =1 

HFI (1) “NBiADE 

ALX (2) =AUX (1 ) 

AUX(2)=PAI/4G. DO 

NHN (2) =T 

NST (2) =1 

NFI (2)=NBLADE 

ALX (3) = AUX (2) 

AUX (3) =PAI/4# DO 
NHN (3) =1 
NST (3) =1 
NFI (31 =NBLADE 
ALX (4) =AUX (3) 

AUX (4) =PAI 
NHN (4) =1 
NST (4) =1 
NFI (4) =N BLADE 
ALX (5) =AUX (4) 

AUX (5) =2*D0*PAI 
NHN (5) =1 
NST (5) =1 
NFI (5) =NBLADE 
ALX (6) =AUX(5) 
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FL0=F10AT (L) 

AOX (6) =2. D0*PAI*FLO 
NEN(6)=L-1 
NST (6) =1 
NFI (6)=NBLADE 
C 

C*** ehd of input 

c 

FLO=FLOAT (NB) 

HH=2, DO *PA I/FLO 
DO 20 1=1 >N BLADE 
FLO=FLOAT (1-1) 

ZB=HH*FLO 
CB (I)=DCOS (ZB) 

SB (I) =DSIN(ZB) 

20 CONTINUE 

DO 105 INTEG=IBEGIN, IFINAL 
NSTAET=N ST (INTEG) 

NFINAL=NFI (INTEG) 

NH=NHN (INTEG) 

AL=ALX (INTEG) 

AO=AUX (INTEG) 

IF (ICHK oNE, 1). GO TO 100 
WHITE (6,6G) 

60 FORMAT (5X, "INPUT' DATA: "//) 

WHITE (6,61) AL 

61 FORMAT (5X, "LOWER LIMIT =”,F25„t5/) 

WRITE (6,62) AD 

62 FORMAT (5X,"UPER LIMIT =",F25.15/) 

WRITE (6,63) NB 

63 FORMAT (5X, "NUMBER OF BLADES: " ,14/) 

WRITE (6,64) NSTART, NFINAL 

64 FORMAT (5X, "NSTART =" ,12, 2X, "NFINAL =",I3/) 
WRITE (6,65) NH 

65 FORMAT (5X, "NUMBER OF INTERVALS, NH=»,I3/) 
WRITE (6 ,42) 
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4 2 FORMAT (5X,» BLADE CONSIDERED : »»/) 

WRITE (6,43) (I,I=NSTART, NFINAL) 

4 3 FORMAT (5 X, " BLADE ”,12/) 

IF (ICHK.EQ. 1) WRITE (6,26) 

26 FORMAT (5X,"NH«,5X,»R",2X,"ROU" ,2 X, "RAMDAK"/) 

100 CALL GQDDAD (AL,AU) 

DO 70 K=1,NPRO 
ROU=RQ (K) 

RAMDAK=RAKDA (K) 

DO 71 J=1,NPR 
SR=R { J) 

ANS=0. DO 

IF (DABS (RR-ROO) • LE. t. D-10) GO TO 72 
ANS=GQUZ25 (FAXIAL) 

7 2 AXIAL (J , K) =AXIAL (J ,K) + ANS 

IF (ICHK.EQ, 1) WRITE (6,15) NH,RR ,ROU ,RAHDAK ,ANS 
7T CONTINUE 
70 CONTINUE 

15 FORMAT (5X,I2,5X,F10, 7,5X ,F10. 7 , 2X, FIG. 7, 2X, D35, 2 8/) 
105 CONTINUE 
RETURN 
END 
C 
C 

FUNCTION GQUZ25 (AUX) 

C 

C PURPOSE 

C 

C TO COMPUTE IA12 FOR A GIVEN INTERVAL BY USING 

C 25-POINT GAUSS-LEG EN DR E FORMULA 

C 

IMPLICIT DOUBLE PRECISION (A-H ,0-Z) 

DOUBLE PRECISION W (25) 

COMMON /GQUZ5 /NH ,N START, NFINAL 
COMMON /GQUZ6/HLETH 

COMMON /GQUZ7/ZZ (25,5) ,CZZ (25 ,5) ,S ZZ (25 ,5) 



W (1 > =0,1 13937985010262 87 S479029641 132D-1 
W (2) =0. 2635 498661503 2137261901 8152953D-1 
W (3) =0 • 40939156701 30 63 12 65 56 23487711 6D-1 
W (4)=C. 5490 4695975 8351 91 925936891 5405D-1 
H (5) =0* 6803 83338123569 1 72071 871856567D-1 
W (6} =G. 8014 07003350C 10 1801 3234959 669 1 D- 1 
W (7) = 0« 91028261 982 96364981 14972207029D-1 
W (8) =0*. 100 5 3 59 490670 5Q6442O22O6890393DO 
W {9) =0,1085 1962 447 42 63 6531160 93957050D0 
W (10) =C.» 11 48582591 4571 164833 9325545870D0 
W (11) =0.119455763535784772228178126513D0 
W (12) -0» 12224244299031004168895951 8946D0 
W (13) =0. 1 23 17 60 5 3726 71 54 5 120 39 02 87307 9D0 
DO 1 1=14,25 
« (I)=W (26-1) 

CONTINUE 

H=HLETH 

NST=NSTART 

NEI=NFINAL 

SUM=0.D0 

DO 9 1=1,25 

H=0* DO ■ 

DO 10 K=1,NH 
Z=ZZ(I,K) 

CZ=CZZ (I,K) 

SZ=SZZ(I,K) . 

R=R+AUX (NST,NFI,Z,CZ,SZ) 

CONTINDE 
SUM=SUH+W (I) *R 
CONTINDE 

GQOZ25=0.5DO *H*SUM 

RETURN 

END 


SUBROUTINE GQUDAD (A,B) 



PURPOSE 


TO COMPOTE THE COEFFICIENTS OSED IN FUNCTION GQUZ25 

IMPLICIT DOUBLE PRECISION (A-E ,0-Z) 

DOUBLE PRECISION T {25) 

COMMON /GQUZ5/NH,N START, NFINAL 
COMMON /GQUZ6/HLBTH 

COMMON /GQUZ7/ZZ (25,5) ,CZZ (25,5) ,SZZ (25,5) 

T(1)=-0 S 9955569 697 90 49 80 97 90 87 84 94 6 89 4D0 
T (2) =-0,97666392145951751149831538648000 
T (3) =-0. 942974571228 97 4339414011 1 6 9658D0 
T (4) =-0* 89 4 99 19 978 78 275368 85 10 4200 67 83 D0 
T (5) =-0„ 833442628760 83 4001 421C21 1G8694D0 
T (6) =-0,7592592630 37 357630 5772 8286 5204 DO 
T (7) =— 0* 6735663684734683644851 20 6 3 3248 D0 
T (8) =-0» 57766 29302 41 2229677236 89841 61 3 DO 
T (9) =-0,47300273144571496052218211500900 
T (1 0) =— 0* 36 1172 305 80 93 87 83 7735 82 1730 128 DO 
I (1 1) — 243866 883720 9 88 43 2045 19 03627 97 DO 
T (12) =-0,12286 469261071039638735981880800 
T (13) =0, DO 
DO 1 1=14,25 
T (I) =-T (26-1) 

CONTINUE 

H= (B-A) /FLOAT (NH) 

HLETH=H 

FORMAT (2X,D31, 24,22,12) 

DO 2 1=1,25 
DO 3 K=1,NH 

Z=A+.5*FLOAT (2*K-1) *H+, 5*H*T (I) 

ZZ (I,K)=Z 

CZZ (I , K) =BCOS (Z) 

SSZ (I,K)=DSIN (Z) 
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3 CONTINUE 
2 CONTINUE 
RETURN 
END 
C 
C 

FUNCTION FAXIAL (NST,NFI, Z,CZ,SZ) 

c 

C PURPOSE 

C 

C TO COMPUTE THE INTEGRAND OF IA12 

C 

IMPLICIT DOUBLE PRECISION (A-H ,0-Z) 

COMMON /BLADE/CB (10) ,SB(1G) 

COMMON /AD? AN/ BAM DAK 
COMMON /COOBD/RR ,ROU 
SUM=0» DO 
DO 1 K=NST, NFI, 

C=CZ*CB (K)-SZ*SB (K) 

A=- (ROU-RR) *ROU* (ROU-RR*C) 
B=ROU*ROU+RR*RR-2, DQ*ROU*RR*C+RAMD AK* 

$ RAMDAK*Z*Z 
SUM=SUM+ A/DSQST (B*B*B) 

1 CONTINUE 
FAXIA L=SUH 
RETURN 
END 
C 
C 

SUBROUTINE INDA3 (EPS ,N PK , N PRO , R , RO , RA MDA , 
$ AXIAL, KR,KRO) 

C 

C PURPOSE 

C 
C 
C 


TO COMPUTE IA3 
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C 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
REAL EPS f R (KS) , HO (KSO) , R AM DA (KR,KRO) 
REAL AXIAL (KR,.KRO) 

COMMON /COO RD/RR,ROU 
COMMON /ADVAN/RAHDAK 
EPSLOR=EPS 
DO 1 K=1,NPRO 
ROU=RO (K) 

RAMDAK^RAMDA (K) 

DO 2 J=1 t NPR 
RR=R(J) 

ANS=0o DQ 

IF (DABS (RB- ROD) « LE. 1»D-10) GO 10 3 
ANS=AIND3 (EPSLON) 

3 AXIAL (J,K)=AXIAL (J,K) +ANS 
2 CON TIN DE 
1 CONTINUE 
RETURN 
END 
C 
C 

FUNCTION AIND3 (EPSLON) 

C 

C PURPOSE 

C 

C TO COMPUTE IA3 

C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON /COO BD/RR # ROU 
COMMON /AD VAN /RAM DAK 
A= (ROU-RR)* (ROU-RR) 

B= RR*ROU+RAMDAK*RAMDAK 
C=DSQRT (A+B*EPSLON*EPSLON) 

SUM— BOU*EPSLQN /C 
AIND3=SUM 
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IF (DABS (A) , LB* 1 « D-15) RETURN 
F=G* 5D0*ROU*RR* (SOU- RE) 

U=-EPSLON/(B*C) 

V=DLOG ( (EPSLOH*DSQST (B) + C) /DABS (ROU-RR) ) 

SUM=SUM+F* (U+7/DSQRT (B*B*B) ) 

AIND3=SUH 

RETURN 

END 

C 

C 

SUBROUTINE INDA4 (NPR,NPRO, R,RO ,RAHDA, AXIAL, KR,KSO) 
C 

c 

C PURPOSE 

c 

C TO COMPUTE IA4 

C 

C 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

REAL R (KR) , BO (KRO) , RAM DA (KRO) , AXIAL (KR , KRO) 

COMMON /IND4/AN (15) ,C (15,15) , COS3K (1 5,1 4) ,COS31 
COMMON /COO BD /R R , SO U 
COMMON /ADVAN/RAHDAK 
REWIND 1 

READ (1 ,1)NBLADE,NTERMS,LUP3R 

1 FORMAT (3 (2X,I3) ) 

PAI=3 V 141 59 2653 5 89 79 32 38 46 2643383 27 9 DO 
READ (1 ,2}CQS31 . . ■ 

2 FORMAT (2X,D21, 1 4} 

DO 3 I-1,NTEBH5 
KJ=I+1 

READ (1,4) (COS3K (I, J) ,J-1,KJ) 

3 CONTINUE 

4 FORMAT (5 (2X,D21. 14)) 

DO 5 1=1 , NTERMS 

READ (1,6) (C (I,J) ,J=1 ,1) 
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5 CONTINUE 

6 FORMAT (3 (2X,D35« 28) ) 

L=LUP ER 
FLO=FLOAT (L) 

UPER=2. D0*PAI*FLO 
NB=NBLADE 
CHK=100« DO 

DO 7 K=1 , NPEO 
ROU=RO (K) 

RAMDAK=RAMDA (K) 

I? {DABS (RAMDAK-CHK) , GE. 1*D-8) CALL COEAN (UPER,NTERMS) 
DO 8 J=1 , NPR 
RR=R (J) 

ANS=G» DO 

IF (DABS (RR-ROU) • LE» 1 a D— 1 0) GO TO 9 
ANS=AIND4 (UPER, N TERMS, NB) 

9 AXIAL (J,K)=AXIAL (J,K) +ANS 
8 CONTINUE 

7 CONTINUE 
RETURN 
END 

C 

c 

c 

FUNCTION AIND4 (UPER, NTESMS ,NB) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE IA4 

C 

c - 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

REAL FLOAT 

COMMON /IND4/AN (15) ,C (15,15) ,COS3K (15,14) ,CDS31 
COMMON /COO RD/RR,ROB 
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COMMON /ADVAN/RAHDAK 
U= 0 PE R 

RAMDAO=RAMDAK 

X=RR 

XO=ROU 

P=0. 25 DO* {X *X+X0 *XO) 

Q=~G» 5D0*X*XO 
F=XO*U 
B=FLO AT (NB) 

CHK=0 * DO 

StJM=F* (G*5DO*XO*B-X*COS31) /(2AMDA0*U) **3 
ANS=SUM* (20-X) 

AIND4=ANS 

IF (NTERf!S»LE*0 ) RE TORN 

NMAX=NTERMS 

DO 1 N=1,MAX 

T=F*P**{N) * (XO*B/FLOAT(2*N+2)-X*COS3K (N,1) ) 
KMAX=N 

DO 2 K— 1 , KMAX 

T=T+C (N,K)*P**(N-K)*Q** (K) *F*(XO*COS3K (N , K) - 
$ X*COS3K (N,K+1) ) 

2 CONTINUE 
T=T*AN (N) 

SUM=SUtl+T 

IF (DABS (SUM-CHK) . LE. 1, D-16) SO TO 3 
CHK=SUM 
1 CONTINDE 

3 AIND4=SUM* (XO-X) 

RETURN 

END 

C 

C 

C 

SUBROUTINE COEAN (UPEK,NTERM5) 

C 

C 
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C PURPOSE 

C 

C TO COMPUTE AN (N) 

C 

C 

IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

COMMON /IND4/AN (15) ,C{15,15) ,COS3K (15,14) ,C0S31 

COMMON /ADVAN/RAMDAK 

RAMDAO=RAMDAK 

U— UPER 

NMAX^NTERMS 

A=DLOG (RAMDAO*U) 

SIGN=-1» . 

R=2.*DLOG (2. DO) 

S=DLGAMA (1,1) 

DO 1 N=1,NMAX 

T= DLG AHA (N+ 1,.1) +FLOAT(N) *R- (FLOAT (2*N+3) * 

$ A+S+DLGAMA (N+1 ,0) ) 

AN (N)=SIGN*DEXP (T) 

SIGN=-SIGN 
1 CONTINUE 
RETURN 
END 
C 
C 
C 

FUNCTION DLGAMA (M r N) 

C 

C 

C PURPOSE 

C 

C 

C TO CALCULATE THE NATURAL LOGARITHM OF 

C GAMMA FUNCTION 

C 
C 
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C 

C 

C 

C 

c 

c 

c 

c 


DESCRIPTION OF PARAMETERS 

M INTEGER 

N CONTROL PARAMETER 

IF N=1, DLGAMA=LOG (GAMMA (M+ 1/2) 

IF N * N Ea " 1 , DLGAMA = LOG (GAMMA (M) ) 


IMPLICIT DOUBLE PRECISION (A-H ,0-Z) 

X= FLOAT (M) 

IF (N ,EQ. 1) GO TO 1 
IF (M . GE, 1) GO TO 2 
PRINT 3,H 

3 ‘ FORMAT (3X,"*** AN ERROR OCCURS IN DLGAMA 
$ DUE TO M='» jr I10 ,»***») 

STOP 

2 DLGAMA=O.DQ 

IF ( M* EQ* 1 ) RETURN 
SUM=D*D0 „ 

N MA X= H.~ 1 
Z=G.D0 

DO 4 I=1,NMAX 
Z=Z+1,D0 
SUM=SUM+DLOG (Z) 

4 CONTINUE 
DLGAMA=SUK 
RETURN 

1 NHAX=M 

H=0*5D0 *1* 144729885849400 1741 43427351 353D0 

SUM=0» DO 

Z=0.DQ 

DO 5 1=1 , NM AX 
Z=Z+1„ 

SUH=SUH+DLOG (Z-0.5DG) 

5 CONTINUE 
DLGAMA=H+SUM 
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RETURN 
END 
C • 

C 

C 

SUBROUTINE WHEN? (NB , B , RO ,RAMDA K r AXIAL, TA NG) 
C 
C 

C PURPOSE 

c 

C TO COMPUTE THE AXIAL AND TANGENTIAL 

C INDUCTION FACTORS AT POINT (R,RO) BY 

C USING WRENCH MODIFIED FORMULAS 

C 

C DESCRIPTION OF PARAMETERS 

C 

C NB 

C X 

C XO 

C RAMDAO 

C AXI'AL 

C TANG 

C 
C 

IMPLICIT DOUBLE PRECISION (A-H,0“Z) 

REAL S, RO, RAM DAK, TANG, AXIAL 
X=R 
XO=RO 


NUMBER OF PROPELLER BLADES 
RADIAL COORDINATE OF THE REFERENCE POINT 
RADIAL COORDINATE OF THE SOURCE POINT 
HYDRODYNAHICAL ADVANCE COEFFICIENT AT XO 
AXIAL INDUCTION FACTOR 
TANGENIAL INDUCTION FACTOR 


EAMDAO=RAHDAK 

IF (DABS (X“XO) • GT* 1 • D-IO) GO TO 2 

U=DSQRT (RAHDAO*RAHDAO+X*XO) 

AXIAL=XO/U 

TANG=-RAMDAO/U 

RETURN 

2 Y=X/RAMDAO 
YO=XO/RAMDAO 
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A=DSQET (I* D0+Y**2) 

B=DSQET(1. D0+Y0**2) 

C— (YO* ( A~1 , DO) / (Y* (B-1.D0))) ** {SB) 

D= (A-B) AFLOAT (NB) 

U=C*DEXP(D} 

A-DSQET (DSQR3 ( {1 , 00+ YO **2) / (1 * DO +Y**2 ) ) ) / 

$ (FLOAT (2*NB)*Y0) 

E= (9. D0*YO**2+2,DQ)/DSQSI( (1.DO+YO**2)**3) 
C* (3* D0*Y**2-2i 00) /DSQET ((1, D0+Y**2) **3) 

D= (B+C) /FLOAT (24*NB) 

E*1. DO / (1* DG/U-1, DO) 

F=1.D0/(U-1.D0) 

IF (Y . LT. . 10) F1=-A* (E+D*DLOG (1, DO + E) ) 

IF(Y * GT* YO) F2— A* { F-D*DL03 (1 #DO+F) ) 

IF (Y, GT. YO) GO TO 1 

FLO=FLOAT (NB) 

AXIAL = F10*Y0* (1.DO-Y/IO)* (1,00-2* DO *FLO* 
$ Y0*F1) 

TANG=-2,D0*FLO*FLO*YO* (1 ,DO-YO/Y) *F1 
RETURN 

1 FLO=FLOAT (NB) 

AXIAL=2* DC*FLO*FLO*IO*Y* (1 . DO-YO/Y) *F2 
TANG=-FLO* (1* DO- YO/Y) * (1.00+2* D0*FLO* Y0*F2) 
RETURN 
END 
C 

C ************************************* 

c ***** ***** 

C ***** ACOUSTIC MODEL ***** 

C ***** ***** 

C ************************************* 

C 

c 

SUBROUTINE MEAN SQ 
C 

C PURPOSE 
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C 

C • TO COMPUTE- THE ENSEMBLE MEAN SQUARE OF 
C THE PROPELLER ACOUSTIC NOISE DUE TO 

C ROTATIONAL FORCES 

C 

C HEIGHTS OF A 25- POINT GAUSS- LEGENDRE 

C FORMULA 

C 

COMMON /DATA24/SQMEAN 
COMMON /ACOUST/PRESU {25) 

DOUBLE PRECISION W (1 3) 

W (1) =0.1 139379850 1026287 94790 29641 132D-1 
W(2) =0.2635498661503213726190181 52953 D-1 
H (3) =0, 40939156701 30 631265562348771 16D-1 
W (4) =0. 54904695 9758351919259358915405D-1 
W (5) =0. 6 803833381235691 72071 871 856567D-1 
W (6) =0.8014 070033500 10 1801323 49596691 D-1 
H (7) =0.9102 8261 982 96 36498114972207029D-1 
W (8) =0# 1005359 490670 50 64 42022068903 93 DO ' 
W (9) =0.1085 1 962447 426365311 60 93957050D0 
W (10) =0,1 148582591 4571 1 64833 932554587GD0 
W (11) =0. 119455763535784772228178126513D0 
H (12) =0.12224244 2990310041 68895951 8946 DO 
H (13) =0. 12317605372671 5451 20 39 02 873O79D0 
CALL PRESUF 
SUM=0. . 

DO 1 1=1,25 
Z=PRESU (I) 

IF (I. LB. 13) SUM=SUM+S (I) *Z*Z 
IF (I. GT.13) SUM=SUM+H(26-I)*Z*Z 
1 CONTINUE 

SQMEAN=SUM*0.5 

RETURN 

END 

C 

c 
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SUBROUTINE PRESUF 
C 

C PURPOSE 

C 

C TO COMPUTE’ THE PROPELLER' ACOUSTIC 

C PRESSURE 

C 
C 

COMMON /CON ST/PAI 
COMMON /DATA6/HUB 
COMMON /DAT A1 1/MCIRCU 
COMMON /DATA14/G (15) 

COMMON /DATA20/N NOISE, NO SCHB 

COMMON /DATA 22/P KO (25,31 ) ,PKA (25,31) ,PKT (25,31 
COMMON /DATA31/X VELDG,IVELCHB 
COMMON /DATA33/UACHB (51) ,UTCHB(51) 

COMMON /ACOUST/PRESU (25) 

DO 1 1=1 ,25 
PRESU (I) =0. 

1 CONTINUE 

PACT= ( 1 * “HUB) *PAI/8. • 

DO 2 M=1 ,MCIRCU 
IF ( (H—1) «,GT. NNOISE) GO TO 3 
IF ( (H + 1) . GT,NNOISE) GO TO 4 
DO 5 1=1,25 

PRESU (I) =PRSSU(I) +FACT+G (H) * {PKO (I,M) - 
S PKO (I , M+2) ) 

5 CONTINUE 
GO TO 2 

4 DO 6' 1=1,25 

PRESU (I) =PRESU (I) +FACT*G (M) *PKG (I,M) 

6 CONTINUE 

2 CONTINUE 

3 FACT* (1. -HUB) *PAI/16, 

DO 8 1=1,25 

PRESU (I) =PRESU (I)+0„ 25*FACT*G (1) * (PKA (1,1) * 
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$ UACHB (1) +PKT(I,1)*UTCHB (1 ) ) 

8 CONTINUE 
WA=FACT*UACHB (1)/2. . 

WT=FACT*UTCHB (1)/2, 

DO 9 H=1,MCIRCU 

IF (M* NE« 1 ) GO TO 10 
DO 11 1=1,25 

PRESU (I) =P3ESU (I) + G(M) * (HA* (PKA (I,M)/2,- 
$ PKA (I,H+2) ) + WT* (PKT (I ,H)/2* -PKT (I,M+2) ) ) 

11 CONTINUE 
GO TO 9 

10 IF ( (H-1 ) * GT« NNOISE) GO TO 12 
IF ( (H + 1) . GT. ’ NNOISE) GO TO 13 
DO 14 1=1,25 

PRESU (I)=PHESU (I) + G (H) * (HA* (PKA (I,H) - 
$ PKA (I , M+2) ) +WT* (PKT (I,H) -PKT (1,8+2) } ) 

1 4 CONTINUE 
GO TO 9 

13 DO 15 1=1,25 

PRESU (I) =PRESU (I) + G (H) * (HA*PKA (I,M) +WT*PKT (1,1!) ) 

15 CONTINUE 

9 CONTINUE 

12 KEND=NNOISE+1 

DO 16 H=1,HCIRCU 
DO 17 KX=1 , KEND 
K=KX-1 
HK=1, ' 

IF (K* EQ» 0) HK=2» . 

IF ( (K+H-1 ) a GT* IVELDG) GO TO 24 
HAT=1 • 

IF ( (K+H— 1 ) g EQ» 0 ) HAT=2. 

DO 18 1=1,25 

PRESU (I) =PKESU (I) +FACT*G (M) * (UACHB (K+H) * 

$ PKA (I ,K + 1 ) +UTCHB (K+H) *PKT (I,K + 1) ) / (HK*HAT) 

1 8 CONTINUE 
2 4 IA=IABS (K-H+1) 
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IF (IA* GT, IVELDG) GO TO 1 9 
H£T=1 * 

IF(ia.EQ-G) HAT=2>. 

DO 20 1=1 ,2 5 

PRESU (I) =PRESU (I) +FACT*G (M) * (UACHB (IA+1) * 

$ PKA{I,K+1) +UTCHB (IA+1 )*PKT (I,K+1) )/(HAT*HK) 

20 CONTINUE 

19 IF ( (M + K + 1) . GT. IVELDG) GO TO 21 
DO 22 1=1,25 

PRESU (I) =PRESU (I) -FACT*G (M) * (UACHB (M+K+2) * 

$ PKA (I,K+1) +UTCHB (M+K+2) *PKT (I,K+1) ) /HK 

22 CONTINUE 

21 IA=IABS (K-M-1) 

IF (IA, GT, I7ELDG) GO To' 17 
HAT=1, . 

IF(IA.EQ.O) HAT=2, 

DO 23 1=1,25 

PRESU (I) =PSESU (I) -FACT*G (M) * { UACHB (I £+1 ) 

$ *PK£ (I,K+1) +UTCHB (IA+1) *PKT (I,K+1 ) ) / (HAT*HK) 

23 CONTINUE 
17 CONTINUE 
16 CONTINUE 

RETURN 

END 

C 

C 

SUBROUTINE PKOAT 
C 
C 

C PURPOSE 

C 

C TO COMPUTE THE COEFFICIENTS OF THE 

C CHEBYSHEV EXPANSIONS OF KO , K A , AND KT 

C 
C 


COMMON /DAT A6/HUB 
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COMMON /DATA1 9/QNOISE (51 ) 

COMMON /DAT A2Q/NNOIS E, NOSCHB 
COMMON /DATA21/ETA (25) 

COMMON /DAT A22/PK0 (25,31) ,PKA (25, 31) , PKT (25,31) 
COMMON /CONST/PAI 

COMMON /DATA8/NBLA DE, KWENB,LO PIND 
DIMENSION AO (51) ,AA (51 ) , AT (51) ,RD (51) 

NRD=51 

JX=N0SCHB+1 • 

NB=NBLADE 
XB=FLQ AT (NB) 

BB=2.*PAI/XB 

HI=PAI/XB 

C COMPOTE PRO (I,J) , PKA (I, J) , AND PKT(I,J), 

C 1=1(1)25, J=1 (1) NNOISE+1 

DO 1 1=1,25 

C COORDINATE TRANSFORMATION 

PHIO=HI* (1, +ETA.(I) ) 

C COMPOTE 'THE FONCTIONAL VALOES OF PKG, 

C PKA, AND PKT 

DO 2 J=1,JX 
AO (J) =0. 

AA (J) =0. 

AT (J)=0. 

2 CONTINOE 

DO 3 K=1 , NBLADE 
DE1TAK=HB*FL0AT (K-1) 

C CALL SUBROOTJNE RETARD TC GET THE 

C RETARDED DISTANCE 

C 

CALL RETARD (PHI0,D2LTAK,ED, NED) 

DO 4 J=1,JX 
EETA=RD (J) 

ROU=0.5* ( (1 ,-HOB) *QN0ISE (J) + 1.+30B) 

CALL FKOAT (BETA , PHIO ,DELTAK,R00, ANSO , ANSA, ANST) 
AO (J) =A0 (J) +AN50 
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aa (j) «aa(J) +ansa 
AT(J) *AT (J) +ANST 

4 CONTINUE 
3 CONTINUE 

C COMPUTE THE COEFFICIENTS OF KQ , KA, AND KT 

MCHEBI=NOSCHB+1 
NDEG=NNOISE 
NA=51 ■ 

NLETH=51 

JHAZ-NNOISE+1 

CALL CHEBCF (AO ,aCH2BY,SDEG,HLBTH,BO ,BD, NED) 

DO 5 J=1 , UMAX 

IF(J. EQ. 1} PKQ(I,J}=R0 

IF(J. NE,1) PRO {I,J)=RD (J-1) 

5 CONTINUE 

CAIL CHEBCF (AA,MCHEBI , NDEG,NLETH, SO, RD, NED) 

DO 6 J=1,JMAX 

IF { J» EQo 1 ) PXA (I ,J) = R0 

IF (J* N E* 1} PKA (I,J) =RD (J-1) 

6 CONTINUE 

CALL CHEBCF (AT, MCHEBX , NDEG ,NLETH, HO , RD, NED) 

DO 7 J=1 , JM AX 

IF {J# EQ» 1 ) PKT(I,J)=B0 

IF (J» NE, 1 ) PKT (I , J) = HD (J-1 ) 

7 CONTINUE 
1 CONTINUE 

REWIND2 
DO 8 1=1,25 
DO 9 J=1 , UMAX 

WRITE (2, 10) PRO (I, J) , PKA (I, J) ,PKT (I, J) 

9 CONTINUE' 

8 CONTINUE 

10 FORMAT (3 (2X , F23. 16 ) ) 

20 FORMAT (//, 5X, tt*****************************”//) 

21 FORMAT (5X,» REPLACE F.ILS PKO ATD !?»/) 

22 FORMAT (5X , h*****************************”//) 
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WRITE (6,20) 

WRITE (6,21) 

WRITE (6,22) 

RETURN 

END 

C 

c 

SUBROUTINE FKOAT (BETA, PHIO ,DELTAK,ROU, 

$ ANSO ,ANSA, ANST) ■ 

C 

C PURPOSE 

C 

C TO COMPUTE KQ ,KA ,, AND KT 

C 

c 

COMMON /DATA8/NBLADE,KWENH,LUPIND - 
COMMON /DATA12/FM, TM,VF, EAMDAP 

COMMON /DATA17/XDISTA, AZIMUT, ANGLE, SA NGLE, CANGLE 
SI=SIN (PHIO+DELIAK+T H* BETA) 

CI=COS (PHIO+DELTAK+TM*RETA) 

Cl =ROU/RETA- (XDISTA/RETA) *SANGLE*CI 

C2= (XDISTA/RETA) *CANGLE —PM 

C3= (XDISTA/RETA) *SANGLE* SI 

RH=-FM*C2+BOU*TM*C3 

DRM=ROU*TM*Cl 

S QM= FM * FH+ ( RO U* TH* RO U*TM) 

CPLDS=1»-BM 

C4 = (1 . -SQM) /RETA +TM* DR M 
C5=0* 5/ (XDISTA*CPLUS *C PLUS) 

ANSG=C5* (TM*VF*C1 + (C4/CPLUS) * (ROU*C2/ 

$ RAMD AP+7F*C3) ) 

ANSA=C5* (— X DISTANT M* SA NGLE*CI/RETA+ 

$ (C4/CPLUS) *C3) 

ANST— C5* (FM/BETA+ (C4/CPLUS) *C2) 

RETURN' 

END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE RETARD (PHIO, DELTAS, RD,NRD) 


PURPOSE 

TO COMPUTE THE RETARDED DISTANCE 


DOUBLE PRECISION Cl, C2,C3,C4,C5,C6 ,C7 
DOUBLE PRECISION C8,C9 , TRY , FENTON 
DIMENSION RD(NRD) 

COMMON /DATA6/HUB 

COMMON /D AT A8/NB LA DE , K HENH ,1 UPIND 
COMMON /DATA1 2/FM , TM , VF , RAM DAP 

COMMON /DATA17/XDISTA,AZIMUT, ANGLE, SANGLE, CANGLE 
COMMON /DATA1 9/QNOISE (51) 

COMMON /DATA20/NNOISE, NQSCHB 

COMMON /RETIME/C 1, C2 , C3, 04,05,06,07,08, C9 

C1=1,DG-FM*FM 

C2=FM*CANGLE 

C3=C2*C2 

C7=PHI0+DELTAK 

C8=1, DO/XDI ST A 

C9=TM 

IX=N0SCHB+1 

TRY=XDIST A* (-C2+DSQRT {C3+C 1) )/Cl 
DO 1 1=1, IX 

R=0» 5D0* { {1 * DO” HUB) *QNOISE {I}+1* dQ+HUB) 
C4=1.D0+R*R/ (XDISTA*XDISTA) 

C5=2* * R* S AN GLE/X DIST A 
C6=TM*E*SANGLE/XDISTA 
RD (I) =FEWTON (TRY) 


TRY=RD (I) 
1 CONTINUE 
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RETURN 

END 

C 

c 

FUNCTION FEWTON (TRY) 

C 

C 

C PURPOSE 

C 

C TO COMPUTE THE RETARDED DISTANCE BY USING 

C' NEWTON METHOD 

C 
C 

IMPLICIT DOUBLE PRECISION (A-H # 0-Z) 

COMMON /RETIME/Cl,C2,C3,C4,C5,C6,C7,C8,C9 

EPSLON=1.D-16 

RNl=TRY ' 

1 BN=SN1 

C=DCOS (C7+C 9*RN) 

S=DSIN (C7+C9*RN) 

SQT=DSQRT{C3+Cl* (C4-C5*C)) 

U=RN*C8~ (-C2+SQT)/C1 

V=C8-C6*S/SQT 

RN1=RN-U/V 

IF (DABS (RNl-RN) ,GT. EPSLON) GO TO 1 
FENTON=RNl 
RETURN 
END 
C- 
C 

SUBROUTINE ABSIAS 
C 
C 

C PURPOSE 

C 

C TO SUPPLY THE ABSCISSAS OF 
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C A 25-POINT GAUSS-LEGEND1E FORMULA 

C 

DOUBLE’ PRECISION T (13) (1 3) 

COMMON /DATA21/ETA (25) 

W(1) =0.11393798501 0262 879479029641 132D-1 - 
ff (2) =0. 263 5 498661503 21 37261901 81 52953D-1 . 
W (3) =Q. 4093915670130 63 126556 2348771 16D-1 
V (4) =0. 5490469597583 51.91 9259368915405D-1 
W (5) =0. 680383338123569172071871856567D-1 
W (6) =0. 801407003350010 1 80132349596691 D-1 
W (7) =0. 91 02 82619829636 4981 1497220702 9D-1 
W (8 >=0.1 00 5 3 59 490 6 70 50 64 42 02 20 6 8 903 93 DO ' 

S (9) =0. 10851962447 4263 6531 1 6 093 9 57050 DO . 

W (10) =0. 11485825914571 1648 339325545870DO 
N (1 1) =0. 119 45 5763535784772228178126513 DO ’ 
W (12) =0.1.222424429 90 310041 68895951 8946D0 
H(1 3) =0.12317 60 53726 715451 203902 8730 7 9D0 . 
T (1) =-0. 9955569697 90498097908784 946894D0 . 
T (2) =-0. 976663921 4595 17511 498315386480 DO . 
T (3) =-0* 942974571228974339414011 1 69658DQ ! 
T (4) =-0.8949919978782753 68851042006783 DO 
T (5) =-0. 83344262876083 4001 4210 21 108694D0 
T (6) =-G. 7592592630373576305772 8286S204D0 
T (7) =-0. 673566 36 84734683 644851 20 6332 4 8 DO 
T (8) =-0» 577662930241 222967723689841 61 3 DO 
T (9) =-0* 4730027314457149 605221 82115009D0 
T (10) =-0.361172305809387 83 7735821730128DO 
T (11) =-0. 243866 883720988 43 20451 903 6 27 97 DO 
T (12) =-0.1228646 92610710 39 638735981 880 8D0 
T (13) =0* DO 
DO 1 1=1,25 

IF(I.LE.13) ETA (I) =T (I) 

IF (I. GT, 13) ETA (I) =-T (26-1) 

1 CONTINUE 
RETURN 
END 
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PROGRAM COECOS (INPUTyOUTPUT,COSSRA,TAPEU=COSGRA, 
$ TAPE6 “OUTPUT) 


PURPOSE 

TO COMPUTE THE COEFFICIENTS NEEDED IN 
CALCULATING THE LIFTING-LINE INDUCTION 
FACTORS BY MEANS OF THE ALTERNATE METHOD 
DECELOPED IN APPENDIX A 

THESE COEFFICIENTS ARE: 

COS 31 = INTEGRAL OF (G (X ,0) /X**3 ) *DX (FROM 1 TO INFINTY) 

COS3K (N,K) =INTEGRAL OF (G (X,K) /X** (2*N + 3) } *DX 
(FROM 1 TO INFINTY) 

C(I,J) = (I!)/((I-J) ! * (J!) ) 

WHERE 

G (X , K) =SU MMATI ON OF COS (Y (I) ) I=1,2,, e NB 

NB=NUMBER OF PROPELLER BLADES 
Y (I) =X*U+ (2* PA I/NB) * (1-1) 


PRECISION: DOUBLE PRECISION 


DESCRIPTION OF PARAMETERS 

NBLADE NUMBER OF PROPELLER BLADES-DEFINED IN MAIN 
PROGRAM (INPUT) 

NTERMS NUMBER OF TERMS DESIRED IN CALCULATING IA4 
DEFINED IN MAIN PROGRAH (INPUT) 
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C L 

C 

C C 

C 
C 

C ILL OUTPUTS ARB WRITTEN ONTO FILE COSGSA 

C 

C C0SG3A CONTAINS : 

C 

c N BLADE 

C NTERMS 

C COS31 

C COS3K (N ,K) 

C C(N,I) 

C 

c 

IMPLICIT DOUBLE PRECISION (A-H ,0-Z) ' 
COMMON /BBB/C (15,15) 

COMMON /CCC/COS3K (15,16) ,COS31 
REWIND4 

PAI=3, 141592653 58979323 8462643383279DG 
C 

C*** INPUT DATA 
C 

NBLADE=4 

NTERHS=12 

1=4 

C 

C*** END OF INPUT 
C 

NB=NBLADS 
XL=FLOAT (L) 

U=2.D0*PAI*Xi 
WRITE (6, 2G) 

FORMAT (5X, "INPUT DATA"//) 

WRITE (6,21) NB 


U=2 * PAI * L LOWER LIMIT OF 
INTEGRAL IA4 (INPUT) 

C (I,J)= (I!) /{ (I-J) ! * (J!) ) 


20 
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21 FORMAT (-52, “NUMBER OF BLADES: »,I3/} 

WHITE (6,22) 

22 FORMAT (52, ’‘NUMBER OF TERMS USED IN THE POWER SERIES”/) 
WRITE (6,23) NTERMS 

23' FORMAT (52, “NTERMS =",I3/) 

WRITE (6,24) 

2 4 FORMAT (5X,» THE LOWER LIMIT = 2* PAX * L”/) 

WRITE (6,25) L 

25 FORMAT (5X,» WHERE L =«I3//) 

CALL COECNK (NTERMS) 

.CALL COESCF (NB,NTERM S, U}- 
WRITE (4,6) NB,NTERMS, L 

6 FORMAT (3 (2X, 13)) 

CCOS31=COS31 
WRITE (4 ,4) C0S31 

4 FORMAT (2X,-D21* 14) 

DO 1 1—1 , NTERMS 
KJ=I+1 

WRITE (4,5) (COS3K(I,J) ,J=1,KJ) 

1 CONTINUE 

5 FORMAT (5 (22, D21* 14)) 

DO 2 1=1 ,NTERMS 

WRITE (4,7) (C(I,J),J=1,I) 

7 FORMAT (3 (22, D35, 28)) 

2 CONTINUE 
ENDFILE4 
WRITE (6,77) 

77 FORMAT (///) 

WRITE (6,8) 

WRITE (6 ,9) 

WRITE (6,10) 

WRITE (6,9) 

WRITE (6 ,8) 

8 FORMAT (52, »*********************************'’} 

9 FORMAT (52,”***** **♦**”) 

10 FORMAT (5X, ”***** REPLACE COSGRA *****'») 
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STOP 

END 

r 

C 

SUBROUTINE: COESCF{NB,NTEKMS,U) 

c 

C PURPOSE 

c 

C TO COMPUTE COS31 AND CQS3K 

C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /CCC/COS3K (15 ,16) ,COS31 
NMAX=NTERMS 

CALL SCNK (NB,U,1 ,3 ,ANS) 

COS31 =ANS 

DO 1 N=1,NMAX 

KHAX=N+1 

NI=2*N+3 

DO 2 K=1,KMAX 

KI=K 

CALL SCNK{NB,U,KI,NI,ANS) 

COS3K (N,K) =ANS 
2 CONTINUE 
1 ' CONTINUE 
RETURN- 
END 
C 

c 

SUBROUTINE SCNK {NB,U ,KX,HI,ANS)- 
C 

C PURPOSE 

C 

C TO COMPUTE THE FOLLOWING INTEGRAL 

C 

C ANS=INTEGEAL OF (G (X , KI) /X** NI) * DX (FROM 1 

C TO INFINTY) 
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C 

C ' WHERE 

C G(X,KI) = SUMMATION OF (COS (T (I) **KI) , 1=1 , . • , NB 

C ' Y (I) = U*X+2*PAI* (I-D/NB 

C 0= THE LOWER LIMIT OF IA4 

C ’ NB= THE NUMBER OF PROPELLER BLADES 

C 

IMPLICIT DOUBLE PRECISION (A-H, 0- Z) 

REAL FLOAT 
AN S=0 * DO 

IF (HI * LE« 1) RETURN 

KT= 1 

NT=1 

IF ( { (KI/2) *2— KI) * EQ. 0) ST=2 
IF ( ( (NI/2) *2-NI) • EQ# 0 )' NT=2 
K=KI/2 
N=NI/2 

IF (KT, EQ. 1) K=(KI+1)/2 
. IF (NT. EQ, 1 ) N= (NI+ 1) /2 
IF (KT , EQ* .2) GO TO 2 
S=DLG AMA (2*K,0) 

SUM=0« DO 
DO 3 J=1,K 

I? ( ( ( (2*J-1)/NB) *NB- (2*J-1) ) . NE. G) GO TO 3 
A=U*FLOAT (2*J—1) 

CALL SINCOS (N, A , SINE, SINO, COSE, COSO) 

IF (NT .EQ, 1) Z=COSO 
IF (NT »EQ* 2) Z=COSE 

SUH=SUM+Z*DEXP (S-DLGAMA (K-J + 1 ,0) -DLGAMA (K+J,0) ) 

3 .CONTINUE 

FLO=FLOAT (NB) 

ANS= FLO* (SUM/ ( (2, DO) ** (2*K-2) )) 

RETURN 

2 S= DLGAMA (2*K+1 ,0) 

SUM-0, DO 
DO 4 J= 1 , K 
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IP ( ( (2*J/NB) *NB-2*J) , NE.*.0) GO TO 4 

FLO=FLOAT (2*J) 

A= FLO*U 

CALL SINCOS (N,A,SINE,SING, COSE, COSO) 

IF (NT « EQ. 1). Z=COSO 
IF (NT • EQ* 2 ) ' Z=COSE 

SUM=SUH+Z*DEXP(S-DLGAHA (K-J+1,0) -DLGAHA (K+J+1 ,0} ) 
4 CONTINUE 

FLO=FLOAT (NI-1) 

ANS=DEXP(S-2.D0 *DLGAHA (K+1,0) ) *• 5DQ / 

■ $ FLO+SUM 

FLO=FLOAT (NB) 

ANS=ANS*FLO/ ( (2. DO) ** (2*K-1) ) 

RETURN 
! END 
C 
C 

SUBROUTINE SINCOS (N, A, SINE ,SINO , CO SE, COSO) 

C 

C PURPOSE 

C 

C TO CALCULATE THE FOLLOWING INTEGRALS FOR 

C N » GE. 1 AND U . GT, 40 

C 

C SINE (N,A) =INTEGRAL OF (SIN (A*X) /X** (2*N) ) *DX 

C COSE (N,A) ^INTEGRAL OF (COS (A*X)/X** (2*N) )*DX 

C SINO (N,A)=IKTEGRAL OF (SIN (A *X) ** ( 2*N- 1 ) ) * DX 

C COSO (N,A) ^INTEGRAL OF (COS (A*X) ** (2*N-1 ) ) *DX 

C 

IHPLICIT DOUBLE PRECISION (A-H ,0- Z) 

REAL FLOAT 
EPSLON-1.D-15 
SUNR— 1 • DO 
SUHS=1.D0 
SUHT=1^ DO 
CBKR=0«D0 
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CHKS=0, DO 
CHKT=G.DO 
FR=1 * DO 
FS=1» DO 
FT=1 « DO 
KH=C 
KS=0 
KT=0 

X~ FLOAT (2*N“1 ) 

Y=X+1 • 

1 IF(KT.EQ.G) FT=-FT*X*Y/A**2 
' X=Y 

Y=X+1, DO 

IF(KR*EQ,G) FR=-FR*X*Y/A**2 

X=Y 

I=X+1. 

IF(KS.EQ.O)' FS=— FS*X *Y /A **2 
IF(KR„EQ* 0) SUKR=SUMR+FR 
IF (RS « EQs 0} SUHS=SUHS+FS 
XF{KT.EQ,0) SUHT=SUHT+FT 

I? (DABS (StJHT-CHKT) /DABS (SDMT) ,LE, EPSLON) KT=1 

IF (DABS {SOHR-CHKR) /DABS (SONR) a LE» EPSLON) KR=1 

IF (DABS (ST3MS-CHKS) /DABS (SDHS) ,IE, EPSLON) KS=1 

IF ( (KR+KS+KT-3) * EQ. . 0) GO TO 2 

CHKR=SUMR 

CHKS=SUHS 

CHKT=SUHT 

GO TO 1 

2 R=SUHR/'A 
FL=FLOAT (2*N) 

S=StmS* FL/A**2 
T^SUHT/A 

FL= FLOAT (2*N-1) 

U=SUHR* FL/A**2 
SI=DSIN (A) 

CO=DCOS (A) 
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SINE=R*CO+S*SI 
SXNO=T*CO+U*SI 
COSE— S*CO-R*SI 
COSO=U*CO-T*SI 
RETURN 
END 
C 
C 

SUBROUTINE COECNK (NMAX) 

C 

C PURPOSE 

C 

C TO COMPUTE C(N,K) 

C WHERE C (N,K) =N!/( (N-K) ! * K!) 

C 

IMPLICIT DOUBLE PRECISION (A-H ,0- Z) 

COMMON /BBB/C (15,15) 

IF {NMAX. EQ, 0) RETU RN 
DO 1 N=1,NMAX 
A=DLGAMA (N+1,0) 

KMAX=N 

DO 2 K=1,KMAX 

C (N,K) =DEXP (A- DLGAMA (K+1 ,0) -DLGAMA (N-K+1 , 0 ) ) 
2 CONTINUE 

1 CONTINUE 

RETURN 
END 
C 
C 

FUNCTION DLGAMA (M,N) 

C 

C PURPOSE 

C 

C TO CALCULATE THE NATURAL LOGARITHM OF GAMMA 

C FUNCTION 

C WHERE M IS AN INTEGER, N A CONTROL PARAMETER 
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C 

C DLGAHA=LOG {GAMMA (M+1/2) ) IP N=1 

C DLGAMA=LOG (GAMMA (M)) IF N* HE, 1 

C 

IMPLICIT DOUBLE PRECISION (A-H , O-Z) 

X= FLOAT (M) 

IF (N , EQ. 1) GO TO 1 
IF {M . GE. 1) GO TO 2 
PRINT 3, M 

3 FORMAT {3X, M *** AN ERROR OCCURS IN 
S DLGAMA DUE TO M=» , 110 ,"***») 

STOP 

2 DIGAMA=0» DO 

IF (M. EQ* 1) ' RETURN 
SUM=S* DO 
NMAX=H-1 
Z=0* DO • 

DO 4 I=1,NHAX 
Z=Z+1 » DO 
SUM=SUM+DLOG (Z) 

4 CONTINUE 
DLGAHA=SUM 
RETURN 

C N=1 

1 N,HAX=a; 

C ' H=LOG (SQRT (PAI) ) =LOG (GAMMA (1/2) ) 

H=*5D0 *1*T447298858494001 74143427351 353D0 
SUM=0o DO 
Z=0, DO 

DO 5 I=1,NMAX 
Z= Z* 1 g DO 

SUM=SUM+DLOG (Z-0.5DO) 

5 CONTINUE 
DLGA MA=H+SU M 
RETURN 

END 



